1 //
2 // Rice-K fading generator example.
3 //
4 // This example generates correlated circular complex Gauss random variables
5 // using an approximation to the ideal Doppler filter. The resulting Gauss
6 // random variables are converted to Rice-K random variables using a simple
7 // transformation. The resulting output file plots the filter's power
8 // spectral density, the fading power envelope as a function of time, and the
9 // distribution of Rice-K random variables alongside the theoretical PDF to
10 // demonstrate that the technique is valid.
11 //
12
13 #include <stdio.h>
14 #include <stdlib.h>
15 #include <math.h>
16 #include <complex.h>
17 #include <getopt.h>
18
19 #include "liquid.h"
20
21 #define OUTPUT_FILENAME "ricek_channel_example.m"
22
23 // print usage/help message
usage()24 void usage()
25 {
26 printf("%s [options]\n", __FILE__);
27 printf(" h : print help\n");
28 printf(" n : number of samples\n");
29 printf(" f : maximum doppler frequency (0,0.5)\n");
30 printf(" K : Rice K factor, linear, greater than zero\n");
31 printf(" O : mean power, linear (Rice-K distribution 'omega')\n");
32 }
33
main(int argc,char * argv[])34 int main(int argc, char*argv[])
35 {
36 // options
37 unsigned int h_len = 51; // doppler filter length
38 float fd = 0.1f; // maximum doppler frequency
39 float K = 2.0f; // Rice fading factor
40 float omega = 1.0f; // mean power
41 float theta = 0.0f; // angle of arrival
42 unsigned int num_samples=1200; // number of samples
43
44 int dopt;
45 while ((dopt = getopt(argc,argv,"hn:f:K:O:")) != EOF) {
46 switch (dopt) {
47 case 'h': usage(); return 0;
48 case 'n': num_samples = atoi(optarg); break;
49 case 'f': fd = atof(optarg); break;
50 case 'K': K = atof(optarg); break;
51 case 'O': omega = atof(optarg); break;
52 default:
53 exit(1);
54 }
55 }
56
57 // validate input
58 if (K < 0.0f) {
59 fprintf(stderr,"error: %s, fading factor K must be greater than zero\n", argv[0]);
60 exit(1);
61 } else if (omega < 0.0f) {
62 fprintf(stderr,"error: %s, signal power Omega must be greater than zero\n", argv[0]);
63 exit(1);
64 } else if (fd <= 0.0f || fd >= 0.5f) {
65 fprintf(stderr,"error: %s, Doppler frequency must be in (0,0.5)\n", argv[0]);
66 exit(1);
67 } else if (h_len < 4) {
68 fprintf(stderr,"error: %s, Doppler filter length too small\n", argv[0]);
69 exit(1);
70 } else if (num_samples == 0) {
71 fprintf(stderr,"error: %s, number of samples must be greater than zero\n", argv[0]);
72 exit(1);
73 }
74
75 unsigned int i;
76
77 // allocate array for output samples
78 float complex * y = (float complex*) malloc(num_samples*sizeof(float complex));
79
80 // generate Doppler filter coefficients
81 float h[h_len];
82 liquid_firdes_doppler(h_len, fd, K, theta, h);
83
84 // normalize filter coefficients such that output Gauss random
85 // variables have unity variance
86 float std = 0.0f;
87 for (i=0; i<h_len; i++)
88 std += h[i]*h[i];
89 std = sqrtf(std);
90 for (i=0; i<h_len; i++)
91 h[i] /= std;
92
93 // create Doppler filter from coefficients
94 firfilt_crcf fdoppler = firfilt_crcf_create(h,h_len);
95
96 // generate complex circular Gauss random variables
97 float complex v; // circular Gauss random variable (uncorrelated)
98 float complex x; // circular Gauss random variable (correlated w/ Doppler filter)
99 float s = sqrtf((omega*K)/(K+1.0));
100 float sig = sqrtf(0.5f*omega/(K+1.0));
101 for (i=0; i<num_samples; i++) {
102 // generate complex Gauss random variable
103 crandnf(&v);
104
105 // push through Doppler filter
106 firfilt_crcf_push(fdoppler, v);
107 firfilt_crcf_execute(fdoppler, &x);
108
109 // convert result to random variable with Rice-K distribution
110 y[i] = _Complex_I*( crealf(x)*sig + s ) +
111 ( cimagf(x)*sig );
112 }
113
114 // destroy filter object
115 firfilt_crcf_destroy(fdoppler);
116
117 // export results to file
118 FILE * fid = fopen(OUTPUT_FILENAME,"w");
119 fprintf(fid,"%% %s, auto-generated file\n\n",OUTPUT_FILENAME);
120 fprintf(fid,"clear all;\n");
121 fprintf(fid,"close all;\n");
122 fprintf(fid,"\n");
123 fprintf(fid,"h_len = %u;\n", h_len);
124 fprintf(fid,"num_samples = %u;\n", num_samples);
125
126 // save filter coefficients
127 for (i=0; i<h_len; i++)
128 fprintf(fid,"h(%6u) = %12.4e;\n", i+1, h[i]);
129
130 // save samples
131 for (i=0; i<num_samples; i++)
132 fprintf(fid,"y(%6u) = %12.4e + 1i*%12.4e;\n", i+1, crealf(y[i]), cimagf(y[i]));
133
134 // plot power spectral density of filter
135 fprintf(fid,"nfft = min(1024, 2^(ceil(log2(h_len))+4));\n");
136 fprintf(fid,"f = [0:(nfft-1)]/nfft - 0.5;\n");
137 fprintf(fid,"H = 20*log10(abs(fftshift(fft(h,nfft))));\n");
138 fprintf(fid,"figure;\n");
139 fprintf(fid,"plot(f,H);\n");
140 fprintf(fid,"axis([-0.5 0.5 -80 20]);\n");
141 fprintf(fid,"xlabel('Normalized Frequency [f/F_s]');\n");
142 fprintf(fid,"ylabel('Filter Power Spectral Density [dB]');\n");
143 fprintf(fid,"grid on;\n");
144
145 // plot fading profile
146 fprintf(fid,"figure;\n");
147 fprintf(fid,"t = 0:(num_samples-1);\n");
148 fprintf(fid,"plot(t,20*log10(abs(y)));\n");
149 fprintf(fid,"xlabel('Normalized Time [t F_s]');\n");
150 fprintf(fid,"ylabel('Fading Power Envelope [dB]');\n");
151 fprintf(fid,"axis([0 num_samples -40 10]);\n");
152 fprintf(fid,"grid on;\n");
153
154 // plot distribution
155 fprintf(fid,"[nn xx] = hist(abs(y),15);\n");
156 fprintf(fid,"bin_width = xx(2) - xx(1);\n");
157 fprintf(fid,"ymax = max(abs(y));\n");
158 fprintf(fid,"s = %12.4e;\n", s);
159 fprintf(fid,"sig = %12.4e;\n", sig);
160 fprintf(fid,"yp = 1.1*ymax*[1:500]/500;\n");
161 fprintf(fid,"pdf = (yp/sig^2) .* exp(-(yp.^2+s^2)/(2*sig^2)) .* besseli(0,yp*s/sig^2);\n");
162 fprintf(fid,"figure;\n");
163 fprintf(fid,"plot(yp,pdf,'-', xx,nn/(num_samples*bin_width),'x');\n");
164 fprintf(fid,"xlabel('Fading Magnitude');\n");
165 fprintf(fid,"ylabel('Probability Density');\n");
166 fprintf(fid,"legend('theory','data','location','northeast');\n");
167
168 // close output file
169 fclose(fid);
170 printf("results written to %s\n", OUTPUT_FILENAME);
171
172 // clean up allocated arrays
173 free(y);
174
175 printf("done.\n");
176 return 0;
177 }
178
179