Elaboradar  0.1
 Tutto Classi Namespace File Funzioni Variabili Tipi enumerati (enum) Gruppi
FIR_filter.h
1 #include <fftw3.h>
2 #include <cmath>
3 #include <iostream>
4 
5 using namespace std;
6 
7 class FIR_filter
8 {
9 public:
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  }
60 private:
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 };