1 //
2 // mdct_example.c
3 //
4 // Modified discrete cosine transform (MDCT) example. Demonstrates
5 // the functionality and interface for the mdct and imdct routines.
6 //
7
8 #include <stdio.h>
9 #include <stdlib.h>
10 #include <string.h>
11 #include <math.h>
12
13 #include "liquid.h"
14
15 #define OUTPUT_FILENAME "mdct_example.m"
16
main()17 int main() {
18 // MDCT options
19 unsigned int num_channels=64; // MDCT size
20 unsigned int num_symbols=16;
21
22 // filter options
23 unsigned int h_len=21; // filter length
24 float fc=0.01f; // filter cutoff
25 float As=60.0f; // stop-band attenuation [dB]
26 float mu=0.0f; // timing offset
27
28 // derived values
29 unsigned int i,j;
30 unsigned int num_samples = num_channels*num_symbols;
31 float x[num_samples];
32 float X[num_samples];
33 float y[num_samples];
34 for (i=0; i<num_samples; i++) {
35 x[i] = 0.0f;
36 X[i] = 0.0f;
37 y[i] = 0.0f;
38 }
39
40 FILE*fid = fopen(OUTPUT_FILENAME,"w");
41 fprintf(fid,"%% %s : auto-generated file\n\n", OUTPUT_FILENAME);
42 fprintf(fid,"clear all;\n");
43 fprintf(fid,"close all;\n");
44 fprintf(fid,"num_channels = %u;\n", num_channels);
45 fprintf(fid,"num_symbols = %u;\n", num_symbols);
46 fprintf(fid,"num_samples = num_channels*num_symbols;\n");
47
48 // generate window
49 float w[2*num_channels];
50 for (i=0; i<2*num_channels; i++) {
51 #if 0
52 // raised half sine
53 w[i] = sinf(M_PI/(2.0f*num_channels)*(i+0.5f));
54 #elif 0
55 // shaped pulse
56 float t0 = sinf(M_PI/(2*num_channels)*(i+0.5));
57 w[i] = sinf(M_PI*0.5f*t0*t0);
58 #else
59 w[i] = liquid_kbd(i,2*num_channels,10.0f);
60 #endif
61 }
62
63 // generate basic filter
64 float h[h_len];
65 liquid_firdes_kaiser(h_len,fc,As,mu,h);
66 firfilt_rrrf f = firfilt_rrrf_create(h,h_len);
67
68 // generate random signal (filtered noise)
69 for (i=0; i<num_samples; i++) {
70 // generate filtered noise
71 firfilt_rrrf_push(f, randnf());
72 firfilt_rrrf_execute(f, &x[i]);
73
74 fprintf(fid,"x(%4u) = %12.4e;", i+1, x[i]);
75 }
76 firfilt_rrrf_destroy(f);
77
78 // run analyzer, accounting for input overlap
79 float buffer[2*num_channels];
80 memset(buffer, 0x00, 2*num_channels*sizeof(float));
81 for (i=0; i<num_symbols; i++) {
82 // copy last half of buffer to first half
83 memmove(buffer, &buffer[num_channels], num_channels*sizeof(float));
84
85 // copy input block to last half of buffer
86 memmove(&buffer[num_channels], &x[i*num_channels], num_channels*sizeof(float));
87
88 // run transform
89 mdct(buffer, &X[i*num_channels], w, num_channels);
90
91 //mdct(&x[i*num_channels],&X[i*num_channels],w,num_channels);
92 }
93
94 // run synthesizer, accounting for output overlap
95 for (i=0; i<num_symbols; i++) {
96 // run inverse MDCT
97 imdct(&X[i*num_channels],buffer,w,num_channels);
98
99 // accumulate first half of buffer to output
100 for (j=0; j<num_channels; j++)
101 y[i*num_channels+j] += buffer[j];
102
103 // copy last half of buffer to output (only if we aren't at the last symbol)
104 if (i==num_symbols-1) break;
105 memmove(&y[(i+1)*num_channels], &buffer[num_channels], num_channels*sizeof(float));
106 }
107
108 // print results to file
109 fprintf(fid,"w = zeros(1,2*num_channels);\n");
110 for (i=0; i<2*num_channels; i++)
111 fprintf(fid,"w(%3u) = %12.4e;\n",i+1,w[i]);
112
113 fprintf(fid,"x = zeros(1,num_samples);\n");
114 fprintf(fid,"y = zeros(1,num_samples);\n");
115 for (i=0; i<num_samples; i++) {
116 fprintf(fid,"x(%3u) = %12.4e;\n", i+1, x[i]);
117 fprintf(fid,"y(%3u) = %12.4e;\n", i+1, y[i]);
118 }
119 fprintf(fid,"X = zeros(num_symbols,num_channels);\n");
120 for (i=0; i<num_symbols; i++) {
121 for (j=0; j<num_channels; j++)
122 fprintf(fid,"X(%3u,%3u) = %12.4e;\n", i+1, j+1, X[i*num_channels+j]);
123 }
124
125 // plot spectral response
126 /*
127 fprintf(fid,"figure;\n");
128 fprintf(fid,"f = [0:(num_channels-1)]/(2*num_channels);\n");
129 fprintf(fid,"plot(f,20*log10(abs(X')),'Color',[1 1 1]*0.5);\n");
130 */
131
132 // plot time-domain reconstruction
133 fprintf(fid,"t = 0:(num_samples-1);\n");
134 fprintf(fid,"figure;\n");
135 fprintf(fid,"plot(t,x,t-num_channels,y);\n");
136 fprintf(fid,"xlabel('sample index');\n");
137 fprintf(fid,"ylabel('signal output');\n");
138 fprintf(fid,"legend('original','reconstructed',0);\n");
139
140 fclose(fid);
141 printf("results written to %s\n", OUTPUT_FILENAME);
142
143 printf("done.\n");
144 return 0;
145 }
146
147