1 //
2 // fftfilt_crcf_example.c
3 //
4 // Complex FFT-based finite impulse response filter example. This example
5 // demonstrates the functionality of firfilt by designing a low-order
6 // prototype and using it to filter a noisy signal. The filter coefficients
7 // are real, but the input and output arrays are complex. The filter order
8 // and cutoff frequency are specified at the beginning, and the result is
9 // compared to the regular corresponding firfilt_crcf output.
10 //
11 // SEE ALSO: `firfilt_crcf_example.c`
12 //
13
14 #include <stdio.h>
15 #include <math.h>
16 #include <complex.h>
17
18 #include "liquid.h"
19
20 #define OUTPUT_FILENAME "fftfilt_crcf_example.m"
21
main()22 int main() {
23 // options
24 unsigned int h_len=57; // filter length
25 float fc=0.10f; // cutoff frequency
26 float As=60.0f; // stop-band attenuation
27 unsigned int n=64; // number of samples per block
28 unsigned int num_blocks=6; // total number of blocks
29
30 // derived values
31 unsigned int num_samples = n * num_blocks;
32
33 // design filter
34 float h[h_len];
35 liquid_firdes_kaiser(h_len, fc, As, 0, h);
36
37 // design FFT-based filter and scale to bandwidth
38 fftfilt_crcf q0 = fftfilt_crcf_create(h, h_len, n);
39 fftfilt_crcf_set_scale(q0, 2.0f*fc);
40
41 // design regular FIR filter
42 firfilt_crcf q1 = firfilt_crcf_create(h, h_len);
43 firfilt_crcf_set_scale(q1, 2.0f*fc);
44
45 unsigned int i;
46
47 // allocate memory for data arrays
48 float complex x[num_samples]; // input
49 float complex y0[num_samples]; // output (fftfilt)
50 float complex y1[num_samples]; // output (firfilt)
51
52 // generate input signal (noise)
53 for (i=0; i<num_samples; i++)
54 x[i] = randnf() + _Complex_I*randnf();
55
56 // run signal through fft-based filter in blocks
57 for (i=0; i<num_blocks; i++)
58 fftfilt_crcf_execute(q0, &x[i*n], &y1[i*n]);
59
60 // run signal through regular filter one sample at a time
61 for (i=0; i<num_samples; i++) {
62 // run filter
63 firfilt_crcf_push(q1, x[i]);
64 firfilt_crcf_execute(q1, &y0[i]);
65 }
66
67 // destroy filter objects
68 fftfilt_crcf_destroy(q0);
69 firfilt_crcf_destroy(q1);
70
71 // compute error norm
72 float rmse = 0.0f;
73 printf(" %6s : %8s : %8s (%8s), %8s : %8s (%8s)\n",
74 "index",
75 "re{fir}", "re{fft}", "re{err}",
76 "im{fir}", "im{fft}", "im{err}");
77 for (i=0; i<num_samples; i++) {
78 float complex e = y0[i] - y1[i];
79 printf(" %6u : %8.5f : %8.5f (%8.5f), %8.5f : %8.5f (%8.5f)\n",
80 i,
81 crealf(y0[i]), crealf(y1[i]), crealf(e),
82 cimagf(y0[i]), cimagf(y1[i]), cimagf(e));
83
84 // accumulate error
85 rmse += crealf( e*conjf(e) );
86 }
87 // normalize RMS error
88 rmse = sqrtf( rmse/(float)num_samples );
89 printf(" rmse : %12.4e\n", rmse);
90
91 //
92 // plot results to output file
93 //
94 FILE * fid = fopen(OUTPUT_FILENAME,"w");
95 fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
96 fprintf(fid,"clear all;\n");
97 fprintf(fid,"close all;\n");
98 fprintf(fid,"\n");
99 fprintf(fid,"h_len=%u;\n", h_len);
100 fprintf(fid,"n=%u;\n",num_samples);
101 fprintf(fid,"x =zeros(1,n);\n");
102 fprintf(fid,"y0=zeros(1,n);\n");
103 fprintf(fid,"y1=zeros(1,n);\n");
104
105 for (i=0; i<num_samples; i++) {
106 fprintf(fid,"x( %4u) = %12.4e + j*%12.4e;\n", i+1, crealf( x[i]), cimagf( x[i]));
107 fprintf(fid,"y0(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(y0[i]), cimagf(y0[i]));
108 fprintf(fid,"y1(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(y1[i]), cimagf(y1[i]));
109 }
110
111 // plot output
112 fprintf(fid,"t = 0:(n-1);\n");
113 fprintf(fid,"figure;\n");
114 fprintf(fid,"subplot(2,1,1);\n");
115 fprintf(fid," plot(t,real(y0),'-','Color',[0 0.5 0.2],'LineWidth',1,...\n");
116 fprintf(fid," t,real(y1),'-','Color',[0 0.2 0.5],'LineWidth',1);\n");
117 fprintf(fid," xlabel('time');\n");
118 fprintf(fid," ylabel('real');\n");
119 fprintf(fid," legend('fftfilt','firfilt','location','northeast');\n");
120 fprintf(fid," grid on;\n");
121 fprintf(fid,"subplot(2,1,2);\n");
122 fprintf(fid," plot(t,imag(y0),'-','Color',[0 0.5 0.2],'LineWidth',1,...\n");
123 fprintf(fid," t,imag(y1),'-','Color',[0 0.2 0.5],'LineWidth',1);\n");
124 fprintf(fid," xlabel('time');\n");
125 fprintf(fid," ylabel('imag');\n");
126 fprintf(fid," legend('fftfilt','firfilt','location','northeast');\n");
127 fprintf(fid," grid on;\n");
128
129 fclose(fid);
130 printf("results written to %s.\n", OUTPUT_FILENAME);
131
132 printf("done.\n");
133 return 0;
134 }
135
136