1 //
2 // resamp2_crcf_interp_example.c
3 //
4 // Halfband interpolator. This example demonstrates the interface to the
5 // interpolating halfband resampler. A low-frequency input sinusoid is
6 // generated and fed into the interpolator one sample at a time,
7 // producing two outputs at each iteration. The results are written to
8 // an output file.
9 //
10 // SEE ALSO: resamp2_crcf_interp_example.c
11 // interp_crcf_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_interp_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 input 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.074f; // input tone frequency
38 unsigned int num_samples = 64; // number of input 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[ num_samples]; // input array
56 float complex y[2*num_samples]; // output array
57
58 // generate input
59 unsigned int w_len = num_samples - 2*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<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 interpolator
79 resamp2_crcf_interp_execute(q, x[i], &y[2*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,"num_samples=%u;\n", num_samples);
95 fprintf(fid,"delay =%u;\n", delay);
96 fprintf(fid,"w_sum =%12.8f;\n", w_sum);
97
98 for (i=0; i<num_samples; i++) {
99 // save results to output file
100 fprintf(fid,"x(%3u) = %12.8f + j*%12.8f;\n", i+1, crealf(x[i ]), cimagf(x[i ]));
101 fprintf(fid,"y(%3u) = %12.8f + j*%12.8f;\n", 2*i+1, crealf(y[2*i+0]), cimagf(y[2*i+0]));
102 fprintf(fid,"y(%3u) = %12.8f + j*%12.8f;\n", 2*i+2, crealf(y[2*i+1]), cimagf(y[2*i+1]));
103 }
104
105 // plot time series
106 fprintf(fid,"tx = 0:( num_samples-1);\n");
107 fprintf(fid,"ty = [0:(2*num_samples-1)]/2 - (delay+1)/2;\n");
108 fprintf(fid,"figure;\n");
109 fprintf(fid,"subplot(2,1,1);\n");
110 fprintf(fid," plot(tx,real(x),'-s','Color',[0.5 0.5 0.5],'MarkerSize',1,...\n");
111 fprintf(fid," ty,real(y),'-s','Color',[0.5 0.0 0.0],'MarkerSize',1);\n");
112 fprintf(fid," legend('original','interpolated','location','northeast');");
113 fprintf(fid," axis([-delay num_samples -1.2 1.2]);\n");
114 fprintf(fid," grid on;\n");
115 fprintf(fid," xlabel('Normalized Time [t/T_s]');\n");
116 fprintf(fid," ylabel('real');\n");
117 fprintf(fid,"subplot(2,1,2);\n");
118 fprintf(fid," plot(tx,imag(x),'-s','Color',[0.5 0.5 0.5],'MarkerSize',1,...\n");
119 fprintf(fid," ty,imag(y),'-s','Color',[0.0 0.5 0.0],'MarkerSize',1);\n");
120 fprintf(fid," legend('original','interpolated','location','northeast');");
121 fprintf(fid," axis([-delay num_samples -1.2 1.2]);\n");
122 fprintf(fid," grid on;\n");
123 fprintf(fid," xlabel('Normalized Time [t/T_s]');\n");
124 fprintf(fid," ylabel('imag');\n");
125
126 // plot spectrum
127 fprintf(fid,"nfft=max([1024, 2^(1+nextpow2(num_samples))]);\n");
128 fprintf(fid,"g = 1/w_sum;\n");
129 fprintf(fid,"X=20*log10(abs(fftshift(fft(x*g, nfft))));\n");
130 fprintf(fid,"Y=20*log10(abs(fftshift(fft(y*g/2,nfft))));\n");
131 fprintf(fid,"f=[0:(nfft-1)]/nfft-0.5;\n");
132 fprintf(fid,"figure;\n");
133 fprintf(fid,"plot( f,X,'LineWidth',1, 'Color',[0.5 0.5 0.5],...\n");
134 fprintf(fid," 2*f,Y,'LineWidth',1.5,'Color',[0.1 0.3 0.5]);\n");
135 fprintf(fid,"grid on;\n");
136 fprintf(fid,"xlabel('Normalized Frequency [f/F_s]');\n");
137 fprintf(fid,"ylabel('PSD [dB]');\n");
138 fprintf(fid,"legend('original','interpolated','location','northeast');");
139 fprintf(fid,"axis([-1 1 -100 10]);\n");
140
141 fclose(fid);
142 printf("results written to %s\n", OUTPUT_FILENAME);
143
144 printf("done.\n");
145 return 0;
146 }
147