12 fftw_plan plan_forw,plan_back;
15 FIR_filter(
unsigned N,
double* ret)
18 in = (
double*) fftw_malloc(
sizeof(
double)*N);
19 out = (fftw_complex*) fftw_malloc(
sizeof(fftw_complex)*(1+floor(0.5*N)));
22 plan_forw = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE);
23 plan_back = fftw_plan_dft_c2r_1d(N, out, re_in, FFTW_MEASURE);
26 void feed(
double* input)
28 copy(input,input+n,in);
36 fftw_destroy_plan(plan_forw);
37 fftw_destroy_plan(plan_back);
45 fftw_execute(plan_forw);
48 fftw_execute(plan_back);
54 for(
unsigned i=0;i<20;i++)
56 cout<<fixed<<in[i]<<
"\t"<<out[i][0]<<
"\t"<<out[i][1]<<
"\t\t"<<re_in[i]/n<<endl;
61 void frequency_filter()
64 hn[0]=hn[20]=1.625807356e-2;
65 hn[1]=hn[19]=2.230852545e-2;
66 hn[2]=hn[18]=2.896372364e-2;
67 hn[3]=hn[17]=3.595993808e-2;
68 hn[4]=hn[16]=4.298744446e-2;
69 hn[5]=hn[15]=4.971005447e-2;
70 hn[6]=hn[14]=5.578764970e-2;
71 hn[7]=hn[13]=6.089991897e-2;
72 hn[8]=hn[12]=6.476934523e-2;
73 hn[9]=hn[11]=6.718151185e-2;
75 fftw_complex* Hz=(fftw_complex*)fftw_malloc(
sizeof(fftw_complex)*(1+floor(0.5*n)));
78 for(
unsigned i=0;i<1+floor(0.5*n);i++)
82 omega=i*M_PI/(1+floor(0.5*n));
83 for(
unsigned j=0;j<20;j++)
85 Hz[i][0]+=hn[j]*cos(-omega*j);
86 Hz[i][1]+=hn[j]*sin(-omega*j);
88 real=out[i][0]*Hz[i][0]-out[i][1]*Hz[i][1];
89 imag=out[i][0]*Hz[i][1]+out[i][1]*Hz[i][0];