Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
FIR_filter.h
1#include <fftw3.h>
2#include <cmath>
3#include <iostream>
4
5using namespace std;
6
7class FIR_filter
8{
9public:
10 double *in,*re_in;//,*re_in2;
11 fftw_complex *out;
12 fftw_plan plan_forw,plan_back;//plan_back2;
13 unsigned n;
14
15 FIR_filter(unsigned N,double* ret)
16 {
17 n=N;
18 in = (double*) fftw_malloc(sizeof(double)*N);
19 out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*(1+floor(0.5*N)));
20 re_in = ret; // punto re_in (che è il risultato) su un puntatore passato al costruttore che deve avere n double allocati
21
22 plan_forw = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE); // sign is not need for real trasforms
23 plan_back = fftw_plan_dft_c2r_1d(N, out, re_in, FFTW_MEASURE);
24 }
25
26 void feed(double* input)
27 {
28 copy(input,input+n,in); // la copia è necessaria, in deve essere costante per plan_forw ed fftw_MEASURE ne distrugge il contenuto
29 // magari esiste un metodo per evitare la copia, ma non lo conosco
30 // si suppone che input abbia n double. Un numero diverso di n richiede l'istanza di un nuovo oggetto FIR
31 // per poter avere i dovuti plan.
32 }
33
34 ~FIR_filter()
35 {
36 fftw_destroy_plan(plan_forw);
37 fftw_destroy_plan(plan_back);
38
39 fftw_free(in);
40 fftw_free(out);
41 }
42
43 void perform()
44 {
45 fftw_execute(plan_forw);
46 frequency_filter();
47 frequency_filter();
48 fftw_execute(plan_back);
49 }
50
51 void dump()
52 {
53 cout<<endl<<endl;
54 for(unsigned i=0;i<20;i++)
55 {
56 cout<<fixed<<in[i]<<"\t"<<out[i][0]<<"\t"<<out[i][1]<<"\t\t"<<re_in[i]/n<<endl;
57 }
58 cout<<endl<<endl;
59 }
60private:
61 void frequency_filter()
62 {
63 double hn[20];
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;
74 hn[10]=6.8001e-2;
75 fftw_complex* Hz=(fftw_complex*)fftw_malloc(sizeof(fftw_complex)*(1+floor(0.5*n)));
76 double omega;
77 double real,imag;
78 for(unsigned i=0;i<1+floor(0.5*n);i++)
79 {
80 Hz[i][0]=0;
81 Hz[i][1]=0;
82 omega=i*M_PI/(1+floor(0.5*n));
83 for(unsigned j=0;j<20;j++)
84 {
85 Hz[i][0]+=hn[j]*cos(-omega*j);
86 Hz[i][1]+=hn[j]*sin(-omega*j);
87 }
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];
90 out[i][0]=real;
91 out[i][1]=imag;
92 }
93 }
94};