#define LEVEL_N	10	
#define INPUT_N 1024  /*(level_N)^2*/
#define K 4
#include<stdio.h>
#include<string.h>
#include<math.h>
#include<stdlib.h>
#include"fwt1d.cpp"
void main()
{	
	char *a;
	int l;
	l = INPUT_N;
	unsigned int sup = 4;            //数列長
	double *s0,*s1,*w1;
	int j = 0;
	char filename[100]="scaling";		//出力ファイル名
	char string[100];
	char string1[100];
	double p[K] = {0.482962913145,0.836516303738,0.224143868042,-0.129409522551}; /*ドベシィの数列 p_k(N=2)*/
	double q[K];	/*ドベシィの数列 q_k(N=2) */
	for(int i = 0 ; i<K; i++){
		q[i] = p[i^1];
	}
	/*入力信号生成（s）*/
	int alt = 1;
	s0 = new double[l];
	for(int i=0; i < 1024; i++){
		s0[i] = sin(2*M_PI*(double)i/102.4);
	}
	FILE *f_s0;
	FILE *f_s1;
	FILE *f_w1;
	sprintf(string,"scaling_%d.dat",j);
	f_s0 = fopen(string,"w");
	for(int i=0;i<l;i++){
		fprintf(f_s0,"%d %lf\n",i*(j+1),s0[i]);
	}
	fclose(f_s0);
	
	strcpy(string,"scaling.dat");
	strcpy(string1,"wavelet.dat");
	f_s1 = fopen(string,"w");
	f_w1 = fopen(string1,"w");

	for(j=1;j <= LEVEL_N; j++){
		printf("level%d\n",j);	
		if(j > 1){
			delete [] s0;
			l = l / 2;	
			s0 = new double[l];
			//一つ前のレベルのスケーリング係数を次の入力にする
			for(int i = 0; i < l ; i++){
				s0[i] = s1[i];
			}
			delete[] s1;
			delete[] w1;
		}
		s1 = new double[l/2];
		w1 = new double[l/2];
		fwt1d(j,s0,l,p,q,sup,s1,w1);	//1次元高速ウェーブレット変換
		
		/*ファイル出力*/
		/*
		sprintf(string,"scaling_%d.dat",j);
		sprintf(string1,"wavelet_%d.dat",j);
		f_s1 = fopen(string,"w");
		f_w1 = fopen(string1,"w");
		*/
		
		for(int i=0;i<l/2;i++){
			fprintf(f_s1,"%d %lf %lf\n",j,i*pow(2,j),s1[i]);
			fprintf(f_w1,"%d %lf %lf\n",j,i*pow(2,j),w1[i]);
		}
	}
	fclose(f_s1);
	fclose(f_w1);
}