1 //
2 // msresamp_crcf_test.c
3 //
4 // Testing the multi-stage arbitrary resampler
5 //
6
7 #include <stdio.h>
8 #include <stdlib.h>
9 #include <complex.h>
10 #include <math.h>
11 #include <getopt.h>
12
13 #include "liquid.h"
14
15 #define OUTPUT_FILENAME "msresamp_crcf_test.m"
16
17 // print usage/help message
usage()18 void usage()
19 {
20 printf("Usage: msresamp_crcf_test [OPTION]\n");
21 printf(" h : print help\n");
22 printf(" r : resampling rate (output/input), default: 0.117\n");
23 printf(" s : stop-band attenuation [dB], default: 60\n");
24 printf(" n : number of input samples, default: 500\n");
25 }
26
main(int argc,char * argv[])27 int main(int argc, char*argv[])
28 {
29 // options
30 float r=1.2f; // resampling rate (output/input)
31 float As=80.0f; // resampling filter stop-band attenuation [dB]
32 unsigned int nx=400; // number of input samples
33 float fc=0.40f; // complex sinusoid frequency
34
35 int dopt;
36 while ((dopt = getopt(argc,argv,"hr:s:n:")) != EOF) {
37 switch (dopt) {
38 case 'h': usage(); return 0;
39 case 'r': r = atof(optarg); break;
40 case 's': As = atof(optarg); break;
41 case 'n': nx = atoi(optarg); break;
42 default:
43 exit(1);
44 }
45 }
46
47 // validate input
48 if (nx == 0) {
49 fprintf(stderr,"error: %s, number of input samples must be greater than zero\n", argv[0]);
50 exit(1);
51 } else if (r <= 0.0f) {
52 fprintf(stderr,"error: %s, resampling rate must be greater than zero\n", argv[0]);
53 exit(1);
54 } else if ( fabsf(log2f(r)) > 10 ) {
55 fprintf(stderr,"error: %s, resampling rate unreasonable\n", argv[0]);
56 exit(1);
57 }
58
59 unsigned int i;
60
61 // derived values
62 unsigned int ny_alloc = (unsigned int) (2*(float)nx * r); // allocation for output
63
64 // allocate memory for arrays
65 float complex x[nx];
66 float complex y[ny_alloc];
67
68 // generate input
69 unsigned int window_len = (3*nx)/4;
70 for (i=0; i<nx; i++)
71 x[i] = i < window_len ? cexpf(_Complex_I*2*M_PI*fc*i) * kaiser(i,window_len,10.0f,0.0f) : 0.0f;
72
73 // create multi-stage arbitrary resampler object
74 msresamp_crcf q = msresamp_crcf_create(r,As);
75 msresamp_crcf_print(q);
76 float delay = msresamp_crcf_get_delay(q);
77
78 // run resampler
79 unsigned int ny;
80 msresamp_crcf_execute(q, x, nx, y, &ny);
81
82 // print basic results
83 printf("input samples : %u\n", nx);
84 printf("output samples : %u\n", ny);
85 printf("delay : %f samples\n", delay);
86
87 // clean up allocated objects
88 msresamp_crcf_destroy(q);
89
90 //
91 // export output
92 //
93 // open/initialize output file
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,"r=%12.8f;\n", r);
100 fprintf(fid,"delay = %f;\n", delay);
101
102 // save input series
103 fprintf(fid,"nx = %u;\n", nx);
104 fprintf(fid,"x = zeros(1,nx);\n");
105 for (i=0; i<nx; i++)
106 fprintf(fid,"x(%6u) = %12.4e + j*%12.4e;\n", i+1, crealf(x[i]), cimagf(x[i]));
107
108 // save output series
109 fprintf(fid,"ny = %u;\n", ny);
110 fprintf(fid,"y = zeros(1,ny);\n");
111 for (i=0; i<ny; i++)
112 fprintf(fid,"y(%6u) = %12.4e + j*%12.4e;\n", i+1, crealf(y[i]), cimagf(y[i]));
113
114 // output results
115 fprintf(fid,"\n\n");
116 fprintf(fid,"%% plot frequency-domain result\n");
117 fprintf(fid,"nfft=1024;\n");
118 fprintf(fid,"%% estimate PSD, normalize by array length\n");
119 fprintf(fid,"X=20*log10(abs(fftshift(fft(x,nfft)/length(x))));\n");
120 fprintf(fid,"Y=20*log10(abs(fftshift(fft(y,nfft)/length(y))));\n");
121 fprintf(fid,"G = max(X);\n");
122 fprintf(fid,"X = X - G;\n");
123 fprintf(fid,"Y = Y - G;\n");
124 fprintf(fid,"f=[0:(nfft-1)]/nfft-0.5;\n");
125 fprintf(fid,"figure;\n");
126 fprintf(fid,"if r>1, fx = f/r; fy = f; %% interpolated\n");
127 fprintf(fid,"else, fx = f; fy = f*r; %% decimated\n");
128 fprintf(fid,"end;\n");
129 fprintf(fid,"plot(fx,X,'Color',[0.5 0.5 0.5],fy,Y,'LineWidth',2);\n");
130 fprintf(fid,"grid on;\n\n");
131 fprintf(fid,"xlabel('normalized frequency');\n");
132 fprintf(fid,"ylabel('PSD [dB]');\n");
133 fprintf(fid,"legend('original','resampled','location','northeast');");
134 fprintf(fid,"axis([-0.5 0.5 -120 10]);\n");
135
136 fprintf(fid,"\n\n");
137 fprintf(fid,"%% plot time-domain result\n");
138 fprintf(fid,"tx=[0:(length(x)-1)];\n");
139 fprintf(fid,"ty=[0:(length(y)-1)]/r-delay;\n");
140 fprintf(fid,"figure;\n");
141 fprintf(fid,"subplot(2,1,1);\n");
142 fprintf(fid," plot(tx,real(x),'-s','Color',[0.5 0.5 0.5],'MarkerSize',1,...\n");
143 fprintf(fid," ty,real(y),'-s','Color',[0.5 0 0], 'MarkerSize',1);\n");
144 fprintf(fid," legend('original','resampled','location','northeast');");
145 fprintf(fid," xlabel('time');\n");
146 fprintf(fid," ylabel('real');\n");
147 fprintf(fid," grid on;\n\n");
148 fprintf(fid,"subplot(2,1,2);\n");
149 fprintf(fid," plot(tx,imag(x),'-s','Color',[0.5 0.5 0.5],'MarkerSize',1,...\n");
150 fprintf(fid," ty,imag(y),'-s','Color',[0 0.5 0], 'MarkerSize',1);\n");
151 fprintf(fid," legend('original','resampled','location','northeast');");
152 fprintf(fid," xlabel('time');\n");
153 fprintf(fid," ylabel('imag');\n");
154 fprintf(fid," grid on;\n\n");
155
156 fclose(fid);
157 printf("results written to %s\n",OUTPUT_FILENAME);
158
159 printf("done.\n");
160 return 0;
161 }
162