1 //
2 // resamp2_crcf_example.c
3 //
4 // This example demonstrates the halfband resampler running as both an
5 // interpolator and a decimator. A narrow-band signal is first
6 // interpolated by a factor of 2, and then decimated. The resulting RMS
7 // error between the final signal and original is computed and printed
8 // to the screen.
9 //
10
11 #include <stdio.h>
12 #include <complex.h>
13 #include <math.h>
14
15 #include "liquid.h"
16
17 #define OUTPUT_FILENAME "resamp2_crcf_example.m"
18
main()19 int main() {
20 unsigned int m=5; // filter semi-length (actual length: 4*m+1)
21 float bw=0.13f; // input signal bandwidth
22 float fc=-0.597f; // input signal carrier frequency (radians/sample)
23 unsigned int num_samples=37; // number of input samples
24 float As=60.0f; // stop-band attenuation [dB]
25
26 unsigned int i;
27
28 // derived values
29 unsigned int n = num_samples + 2*m + 1; // adjusted input sequence length
30
31 // allocate memory for data arrays
32 float complex x[ n];
33 float complex y[2*n];
34 float complex z[ n];
35
36 // generate the baseband signal (filter pulse)
37 float h[num_samples];
38 liquid_firdes_kaiser(num_samples,bw,60.0f,0.0f,h);
39 for (i=0; i<n; i++)
40 x[i] = i < num_samples ? h[i] * cexpf(_Complex_I*fc*i) : 0.0f;
41
42 // create/print the half-band resampler, with a specified
43 // stopband attenuation level
44 resamp2_crcf q = resamp2_crcf_create(m,0,As);
45 resamp2_crcf_print(q);
46
47 // run interpolation stage
48 for (i=0; i<n; i++)
49 resamp2_crcf_interp_execute(q, x[i], &y[2*i]);
50
51 // clear resamp2 object
52 resamp2_crcf_reset(q);
53
54 // execute decimation stage
55 for (i=0; i<n; i++)
56 resamp2_crcf_interp_execute(q, y[2*i], &z[i]);
57
58 // clean up allocated objects
59 resamp2_crcf_destroy(q);
60
61 // compute RMS error
62 float rmse = 0.0f;
63 for (i=2*m; i<n; i++) {
64 float e = cabsf(x[i-2*m] - z[i]);
65 rmse += e*e;
66 }
67 rmse = sqrtf( rmse / (float)(n-2*m) );
68 printf("rms error : %12.4e\n", rmse);
69
70 //
71 // print results to file
72 //
73 FILE*fid = fopen(OUTPUT_FILENAME,"w");
74 fprintf(fid,"%% %s: auto-generated file\n",OUTPUT_FILENAME);
75 fprintf(fid,"clear all;\n");
76 fprintf(fid,"close all;\n\n");
77 fprintf(fid,"bw=%12.8f;\n", bw);
78 fprintf(fid,"n=%u;\n", n);
79
80 // output results
81 for (i=0; i<n; i++)
82 fprintf(fid,"x(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(x[i]), cimagf(x[i]));
83
84 for (i=0; i<2*n; i++)
85 fprintf(fid,"y(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(y[i]), cimagf(y[i]));
86
87 for (i=0; i<n; i++)
88 fprintf(fid,"z(%3u) = %12.4e + j*%12.4e;\n", i+1, crealf(z[i]), cimagf(z[i]));
89
90 // print results
91 fprintf(fid,"\n\n");
92 fprintf(fid,"nfft=1024;\n");
93 fprintf(fid,"f = [0:(nfft-1)]/nfft - 0.5;\n");
94 fprintf(fid,"X = 20*log10(abs(fftshift(fft(x*bw*2,nfft))));\n");
95 fprintf(fid,"Y = 20*log10(abs(fftshift(fft(y*bw, nfft))));\n");
96 fprintf(fid,"Z = 20*log10(abs(fftshift(fft(z*bw*2,nfft))));\n");
97 fprintf(fid,"plot(f,X,f,Y,f,Z);\n");
98 fprintf(fid,"legend('original','up-converted','down-converted',1);\n");
99 fprintf(fid,"grid on;\n");
100 fprintf(fid,"axis([-0.5 0.5 -120 20]);\n");
101
102 fprintf(fid,"\n\n");
103 fprintf(fid,"t0 = 0:[n-1];\n");
104 fprintf(fid,"t1 = 0:[n*2-1];\n");
105 fprintf(fid,"figure;\n");
106 fprintf(fid,"subplot(3,1,1);\n");
107 fprintf(fid," plot(t0,real(x),t0,imag(x));\n");
108 fprintf(fid," legend('I','Q',0);\n");
109 fprintf(fid," axis([0 n -1 1]);\n");
110 fprintf(fid," ylabel('original');\n");
111 fprintf(fid,"subplot(3,1,2);\n");
112 fprintf(fid," plot(t1,real(y),t1,imag(y));\n");
113 fprintf(fid," axis([0 n*2 -1 1]);\n");
114 fprintf(fid," ylabel('interpolated');\n");
115 fprintf(fid,"subplot(3,1,3);\n");
116 fprintf(fid," plot(t0,real(z),t0,imag(z));\n");
117 fprintf(fid," axis([0 n -1 1]);\n");
118 fprintf(fid," ylabel('interp/decim');\n");
119
120
121 fclose(fid);
122 printf("results written to %s\n",OUTPUT_FILENAME);
123
124 printf("done.\n");
125 return 0;
126 }
127