1 //
2 // rresamp_crcf_example.c
3 //
4 // Demonstration of rresamp object whereby an input signal
5 // is resampled at a rational rate Q/P.
6 //
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <complex.h>
11 #include <math.h>
12 #include <getopt.h>
13
14 #include "liquid.h"
15
16 #define OUTPUT_FILENAME "rresamp_crcf_example.m"
17
18 // print usage/help message
usage()19 void usage()
20 {
21 printf("Usage: %s [OPTION]\n", __FILE__);
22 printf("Resample a signal at a rate P/Q\n");
23 printf(" -h : print help\n");
24 printf(" -P <decim> : decimation (output) rate, default: 3\n");
25 printf(" -Q <interp> : interpolation (input) rate, default: 5\n");
26 printf(" -m <len> : filter semi-length (delay), default: 12\n");
27 printf(" -w <bandwidth>: filter bandwidth, default: 0.5f\n");
28 printf(" -s <atten> : filter stop-band attenuation [dB], default: 60\n");
29 }
30
main(int argc,char * argv[])31 int main(int argc, char*argv[])
32 {
33 // options
34 unsigned int P = 3; // output rate (interpolation factor)
35 unsigned int Q = 5; // input rate (decimation factor)
36 unsigned int m = 12; // resampling filter semi-length (filter delay)
37 float bw = 0.5f; // resampling filter bandwidth
38 float As = 60.0f; // resampling filter stop-band attenuation [dB]
39
40 int dopt;
41 while ((dopt = getopt(argc,argv,"hP:Q:m:s:w:")) != EOF) {
42 switch (dopt) {
43 case 'h': usage(); return 0;
44 case 'P': P = atoi(optarg); break;
45 case 'Q': Q = atoi(optarg); break;
46 case 'm': m = atoi(optarg); break;
47 case 'w': bw = atof(optarg); break;
48 case 's': As = atof(optarg); break;
49 default:
50 exit(1);
51 }
52 }
53
54 // validate input
55 if (P == 0 || P > 1000) {
56 fprintf(stderr,"error: %s, input rate P must be in [1,1000]\n", argv[0]);
57 exit(1);
58 } else if (Q == 0 || Q > 1000) {
59 fprintf(stderr,"error: %s, output rate Q must be in [1,1000]\n", argv[0]);
60 exit(1);
61 }
62
63 // create resampler object
64 rresamp_crcf q = rresamp_crcf_create_kaiser(P,Q,m,bw,As);
65 rresamp_crcf_print(q);
66 float rate = rresamp_crcf_get_rate(q);
67
68 // number of sample blocks (limit by large interp/decim rates)
69 unsigned int n = 120e3 / (P > Q ? P : Q);
70
71 // input/output buffers
72 float complex buf_x[Q]; // input
73 float complex buf_y[P]; // output
74
75 // create signal generator (wide-band noise)
76 msourcecf gen = msourcecf_create_default();
77 msourcecf_add_noise(gen, 0.0f, 0.7f * (rate > 1.0 ? 1.0 : rate), 0);
78
79 // create spectral periodogram objects
80 unsigned int nfft = 2400;
81 spgramcf px = spgramcf_create_default(nfft);
82 spgramcf py = spgramcf_create_default(nfft);
83
84 // generate input signal (filtered noise)
85 unsigned int i;
86 for (i=0; i<n; i++) {
87 // write Q input samples to buffer
88 msourcecf_write_samples(gen, buf_x, Q);
89
90 // run resampler and write P output samples
91 rresamp_crcf_execute(q, buf_x, buf_y);
92
93 // write input and output to respective spectral periodogram estimate
94 spgramcf_write(px, buf_x, Q);
95 spgramcf_write(py, buf_y, P);
96 }
97 printf("num samples out : %llu\n", spgramcf_get_num_samples_total(py));
98 printf("num samples in : %llu\n", spgramcf_get_num_samples_total(px));
99
100 // clean up allocated objects
101 rresamp_crcf_destroy(q);
102
103 // compute power spectral density output
104 float X[nfft];
105 float Y[nfft];
106 spgramcf_get_psd(px, X);
107 spgramcf_get_psd(py, Y);
108
109 // export results to file for plotting
110 FILE * fid = fopen(OUTPUT_FILENAME,"w");
111 fprintf(fid,"%% %s: auto-generated file\n",OUTPUT_FILENAME);
112 fprintf(fid,"clear all;\n");
113 fprintf(fid,"close all;\n");
114 fprintf(fid,"P = %u;\n", P);
115 fprintf(fid,"Q = %u;\n", Q);
116 fprintf(fid,"r = P/Q;\n");
117 fprintf(fid,"nfft = %u;\n", nfft);
118 fprintf(fid,"X = zeros(1,nfft);\n");
119 fprintf(fid,"Y = zeros(1,nfft);\n");
120 for (i=0; i<nfft; i++) {
121 fprintf(fid,"X(%3u) = %12.4e;\n", i+1, X[i]);
122 fprintf(fid,"Y(%3u) = %12.4e;\n", i+1, Y[i]);
123 }
124 fprintf(fid,"\n\n");
125 fprintf(fid,"%% plot time-domain result\n");
126 fprintf(fid,"fx=[0:(nfft-1)]/nfft-0.5;\n");
127 fprintf(fid,"fy=fx*r;\n");
128 fprintf(fid,"figure('Color','white','position',[500 500 800 600]);\n");
129 fprintf(fid,"plot(fx,X,'-','LineWidth',2,'Color',[0.5 0.5 0.5],'MarkerSize',1,...\n");
130 fprintf(fid," fy,Y,'-','LineWidth',2,'Color',[0.5 0 0], 'MarkerSize',1);\n");
131 fprintf(fid,"legend('original','resampled','location','northeast');");
132 fprintf(fid,"xlabel('Normalized Frequency [f/F_s]');\n");
133 fprintf(fid,"ylabel('Power Spectral Density [dB]');\n");
134 fprintf(fid,"fmin = min(fx( 1),fy( 1));\n");
135 fprintf(fid,"fmax = max(fx(nfft),fy(nfft));\n");
136 fprintf(fid,"axis([fmin fmax -100 20]);\n");
137 fprintf(fid,"grid on;\n");
138 fclose(fid);
139 printf("results written to %s\n",OUTPUT_FILENAME);
140 return 0;
141 }
142