1 //
2 // rresamp_rrrf_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 <math.h>
11 #include <getopt.h>
12
13 #include "liquid.h"
14
15 #define OUTPUT_FILENAME "rresamp_rrrf_example.m"
16
17 // print usage/help message
usage()18 void usage()
19 {
20 printf("Usage: %s [OPTION]\n", __FILE__);
21 printf("Resample a signal at a rate P/Q\n");
22 printf(" -h : print help\n");
23 printf(" -P <decim> : decimation (output) rate, default: 3\n");
24 printf(" -Q <interp> : interpolation (input) rate, default: 5\n");
25 printf(" -m <len> : filter semi-length (delay), default: 12\n");
26 printf(" -w <bandwidth>: filter bandwidth, default: 0.5f\n");
27 printf(" -s <atten> : filter stop-band attenuation [dB], default: 60\n");
28 }
29
main(int argc,char * argv[])30 int main(int argc, char*argv[])
31 {
32 // options
33 unsigned int P = 3; // output rate (interpolation factor)
34 unsigned int Q = 5; // input rate (decimation factor)
35 unsigned int m = 12; // resampling filter semi-length (filter delay)
36 float bw = 0.5f; // resampling filter bandwidth
37 float As = 60.0f; // resampling filter stop-band attenuation [dB]
38
39 int dopt;
40 while ((dopt = getopt(argc,argv,"hP:Q:m:s:w:")) != EOF) {
41 switch (dopt) {
42 case 'h': usage(); return 0;
43 case 'P': P = atoi(optarg); break;
44 case 'Q': Q = atoi(optarg); break;
45 case 'm': m = atoi(optarg); break;
46 case 'w': bw = atof(optarg); break;
47 case 's': As = atof(optarg); break;
48 default:
49 exit(1);
50 }
51 }
52
53 // validate input
54 if (P == 0 || P > 1000) {
55 fprintf(stderr,"error: %s, input rate P must be in [1,1000]\n", argv[0]);
56 exit(1);
57 } else if (Q == 0 || Q > 1000) {
58 fprintf(stderr,"error: %s, output rate Q must be in [1,1000]\n", argv[0]);
59 exit(1);
60 }
61
62 // create resampler object
63 rresamp_rrrf q = rresamp_rrrf_create_kaiser(P,Q,m,bw,As);
64 rresamp_rrrf_print(q);
65 float rate = rresamp_rrrf_get_rate(q);
66
67 // number of sample blocks (limit by large interp/decim rates)
68 unsigned int n = 120e3 / (P > Q ? P : Q);
69
70 // input/output buffers
71 float buf_x[Q]; // input
72 float buf_y[P]; // output
73
74 // create wide-band noise source with one-sided cut-off frequency
75 iirfilt_rrrf lowpass = iirfilt_rrrf_create_lowpass(15, 0.7f*0.5f*(rate > 1.0 ? 1.0 : rate));
76
77 // create spectral periodogram objects
78 unsigned int nfft = 2400;
79 spgramf px = spgramf_create_default(nfft);
80 spgramf py = spgramf_create_default(nfft);
81
82 // generate input signal (filtered noise)
83 unsigned int i, j;
84 for (i=0; i<n; i++) {
85 // write Q input samples to buffer
86 for (j=0; j<Q; j++)
87 iirfilt_rrrf_execute(lowpass, randnf(), &buf_x[j]);
88
89 // run resampler and write P output samples
90 rresamp_rrrf_execute(q, buf_x, buf_y);
91
92 // write input and output to respective spectral periodogram estimate
93 spgramf_write(px, buf_x, Q);
94 spgramf_write(py, buf_y, P);
95 }
96 printf("num samples out : %llu\n", spgramf_get_num_samples_total(py));
97 printf("num samples in : %llu\n", spgramf_get_num_samples_total(px));
98
99 // clean up allocated objects
100 rresamp_rrrf_destroy(q);
101 iirfilt_rrrf_destroy(lowpass);
102
103 // compute power spectral density output
104 float X[nfft];
105 float Y[nfft];
106 spgramf_get_psd(px, X);
107 spgramf_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