1 //
2 // resamp2_crcf_decim_example.c
3 //
4 // Halfband decimator. This example demonstrates the interface to the
5 // decimating halfband resampler. A low-frequency input sinusoid is
6 // generated and fed into the decimator two samples at a time,
7 // producing one output at each iteration. The results are written to
8 // an output file.
9 //
10 // SEE ALSO: resamp2_crcf_interp_example.c
11 // decim_rrrf_example.c
12 //
13
14 #include <stdio.h>
15 #include <stdlib.h>
16 #include <getopt.h>
17 #include <complex.h>
18 #include <math.h>
19
20 #include "liquid.h"
21
22 #define OUTPUT_FILENAME "resamp2_crcf_decim_example.m"
23
24 // print usage/help message
usage()25 void usage()
26 {
27 printf("%s [options]\n", __FILE__);
28 printf(" h : print help\n");
29 printf(" n : number of output samples, default: 64\n");
30 printf(" m : filter semi-length, default: 5\n");
31 printf(" s : filter stop-band attenuation, default: 60 dB\n");
32 }
33
main(int argc,char * argv[])34 int main(int argc, char*argv[])
35 {
36 unsigned int m=5; // filter semi-length
37 float fc=0.037f; // input tone frequency
38 unsigned int num_samples = 64; // number of output samples
39 float As=60.0f; // stop-band attenuation [dB]
40
41 int dopt;
42 while ((dopt = getopt(argc,argv,"hn:m:s:")) != EOF) {
43 switch (dopt) {
44 case 'h': usage(); return 0;
45 case 'n': num_samples = atoi(optarg); break;
46 case 'm': m = atoi(optarg); break;
47 case 's': As = atof(optarg); break;
48 default:
49 exit(1);
50 }
51 }
52 unsigned int i;
53
54 // allocate arrays
55 float complex x[2*num_samples]; // input array
56 float complex y[ num_samples]; // output array
57
58 // generate input
59 unsigned int w_len = 2*num_samples - 4*m; // window length
60 float beta = 8.0f; // kaiser window factor
61 float w_sum = 0.0f; // gain due to window
62 for (i=0; i<2*num_samples; i++) {
63 // compute windowing function and keep track of gain
64 float w = (i < w_len ? kaiser(i,w_len,beta,0) : 0.0f);
65 w_sum += w;
66
67 // compute windowed complex sinusoid
68 x[i] = w * cexpf(_Complex_I*2*M_PI*fc*i);
69 }
70
71 // create/print the half-band resampler
72 resamp2_crcf q = resamp2_crcf_create(m,0,As);
73 resamp2_crcf_print(q);
74 unsigned int delay = resamp2_crcf_get_delay(q);
75
76 // run the resampler
77 for (i=0; i<num_samples; i++) {
78 // execute the decimator
79 resamp2_crcf_decim_execute(q, &x[2*i], &y[i]);
80
81 // print results to screen
82 printf("y(%3u) = %8.4f + j*%8.4f;\n", i+1, crealf(y[i]), cimagf(y[i]));
83 }
84
85 // destroy half-band resampler
86 resamp2_crcf_destroy(q);
87
88 //
89 // export results
90 //
91 FILE*fid = fopen(OUTPUT_FILENAME,"w");
92 fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
93 fprintf(fid,"clear all;\nclose all;\n\n");
94 fprintf(fid,"h_len=%u;\n", 4*m+1);
95 fprintf(fid,"num_samples=%u;\n", num_samples);
96 fprintf(fid,"delay =%u;\n", delay);
97 fprintf(fid,"w_sum =%12.8f;\n", w_sum);
98
99 for (i=0; i<num_samples; i++) {
100 // save results to output file
101 fprintf(fid,"x(%3u) = %12.8f + j*%12.8f;\n", 2*i+1, crealf(x[2*i+0]), cimagf(x[2*i+0]));
102 fprintf(fid,"x(%3u) = %12.8f + j*%12.8f;\n", 2*i+2, crealf(x[2*i+1]), cimagf(x[2*i+1]));
103 fprintf(fid,"y(%3u) = %12.8f + j*%12.8f;\n", i+1, 0.5f*crealf(y[i ]), 0.5f*cimagf(y[i ]));
104 }
105
106 // plot time series
107 fprintf(fid,"tx = 0:(2*num_samples-1);\n");
108 fprintf(fid,"ty = [0:( num_samples-1)]*2 - delay;\n");
109 fprintf(fid,"figure;\n");
110 fprintf(fid,"subplot(2,1,1);\n");
111 fprintf(fid," plot(tx,real(x),'-s','Color',[0.5 0.5 0.5],'MarkerSize',1,...\n");
112 fprintf(fid," ty,real(y),'-s','Color',[0.5 0.0 0.0],'MarkerSize',1);\n");
113 fprintf(fid," legend('original','decimated','location','northeast');");
114 fprintf(fid," axis([-delay 2*num_samples -1.2 1.2]);\n");
115 fprintf(fid," grid on;\n");
116 fprintf(fid," xlabel('Normalized Time [t/T_s]');\n");
117 fprintf(fid," ylabel('real');\n");
118 fprintf(fid,"subplot(2,1,2);\n");
119 fprintf(fid," plot(tx,imag(x),'-s','Color',[0.5 0.5 0.5],'MarkerSize',1,...\n");
120 fprintf(fid," ty,imag(y),'-s','Color',[0.0 0.5 0.0],'MarkerSize',1);\n");
121 fprintf(fid," legend('original','decimated','location','northeast');");
122 fprintf(fid," axis([-delay 2*num_samples -1.2 1.2]);\n");
123 fprintf(fid," grid on;\n");
124 fprintf(fid," xlabel('Normalized Time [t/T_s]');\n");
125 fprintf(fid," ylabel('imag');\n");
126
127 // plot spectrum
128 fprintf(fid,"nfft=512;\n");
129 fprintf(fid,"g = 1/w_sum;\n");
130 fprintf(fid,"X=20*log10(abs(fftshift(fft( x*g,nfft))));\n");
131 fprintf(fid,"Y=20*log10(abs(fftshift(fft(2*y*g,nfft))));\n");
132 fprintf(fid,"f=[0:(nfft-1)]/nfft-0.5;\n");
133 fprintf(fid,"figure;\n");
134 fprintf(fid,"plot(f, X,'LineWidth',1, 'Color',[0.5 0.5 0.5],...\n");
135 fprintf(fid," f/2,Y,'LineWidth',1.5,'Color',[0.1 0.3 0.5]);\n");
136 fprintf(fid,"grid on;\n");
137 fprintf(fid,"xlabel('Normalized Frequency [f/F_s]');\n");
138 fprintf(fid,"ylabel('PSD [dB]');\n");
139 fprintf(fid,"legend('original','decimated','location','northeast');");
140 fprintf(fid,"axis([-0.5 0.5 -100 10]);\n");
141
142 fclose(fid);
143 printf("results written to %s\n", OUTPUT_FILENAME);
144
145 printf("done.\n");
146 return 0;
147 }
148