1% fdmdv_ut_freq_off.m 2% David Rowe 17 June 2014 3% 4% Unit Test program for freq offset estimation in FDMDV modem. 5% 6% Copyright David Rowe 2012 This program is 7% distributed under the terms of the GNU General Public License 8% Version 2 9 10% [ ] sweep of different delays 11% [ ] sweep of Eb/No 12% [ ] sweep of freq offsets 13% [ ] step change in foff 14% + time to respond 15% [ ] plot/print pass fail/relevant stats 16% + variance 17% + histogram of freq ests? 18 19fdmdv; % load modem code 20hf_sim; % load hf sim code 21 22% --------------------------------------------------------------------- 23% Eb/No calculations. We need to work out Eb/No for each FDM carrier. 24% Total power is sum of power in all FDM carriers. These calcs set the 25% Eb/No of the data carriers, Eb/No of pilot will be higher. 26% --------------------------------------------------------------------- 27 28function [Nsd SNR] = calc_Nsd_from_EbNo(EbNo_dB) 29 global Rs; 30 global Nb; 31 global Nc; 32 global Fs; 33 34 C = 1; % power of each FDM carrier (energy/sample). Total Carrier power should = Nc*C = Nc 35 N = 1; % total noise power (energy/sample) of noise source across entire bandwidth 36 37 % Eb = Carrier power * symbol time / (bits/symbol) 38 % = C *(1/Rs) / Nb 39 Eb_dB = 10*log10(C) - 10*log10(Rs) - 10*log10(Nb); 40 41 No_dBHz = Eb_dB - EbNo_dB; 42 43 % Noise power = Noise spectral density * bandwidth 44 % Noise power = Noise spectral density * Fs/2 for real signals 45 N_dB = No_dBHz + 10*log10(Fs/2); 46 Ngain_dB = N_dB - 10*log10(N); 47 Nsd = 10^(Ngain_dB/20); 48 49 % C/No = Carrier Power/noise spectral density 50 % = power per carrier*number of carriers / noise spectral density 51 CNo_dB = 10*log10(C) + 10*log10(Nc) - No_dBHz; 52 53 % SNR in equivalent 3000 Hz SSB channel, adding extra power for pilot to get 54 % true SNR. 55 56 B = 3000; 57 SNR = CNo_dB - 10*log10(B) + 10*log10((Nc+4)/Nc); 58end 59 60% we keep a m sample buffer in sample_memory 61% update sample_memory with n samples each time this function is called 62% outputs one nfft2 slice of spectrogram in dB. Good idea to make nfft2 a power of 2 63 64function [S, states_out] = spectrogram_update(samples, n, states_in) 65 sample_memory = states_in.sample_memory; 66 m = states_in.m; 67 nfft2 = states_in.nfft2; 68 lower_clip_dB = states_in.lower_clip_dB; 69 dec = states_in.dec; 70 71 sample_memory(1:m-n) = sample_memory(n+1:m); 72 sample_memory(m-n+1:m) = samples; 73 74 F = fft(sample_memory .* hanning(m)', 2*nfft2); 75 S = 20*log10(abs(F(1:dec:nfft2))/(nfft2)); 76 S(find(S < lower_clip_dB)) = lower_clip_dB; % clip lower limit 77 78 states_out = states_in; 79 states_out.sample_memory = sample_memory; 80end 81 82% ------------------------------------------------------------ 83 84function sim_out = freq_off_est_test(sim_in) 85 global Nc; 86 global Nb; 87 global M; 88 global Fs; 89 global pilot_lut_index; 90 global prev_pilot_lut_index; 91 global pilot_lpf1; 92 global Npilotlpf; 93 global spread; 94 global spread_2ms; 95 global hf_gain; 96 97 EbNovec = sim_in.EbNovec; 98 Ndelay = sim_in.delay; 99 frames = sim_in.frames; 100 startup_delay = sim_in.startup_delay; 101 allowable_error = sim_in.allowable_error; 102 foff_hz = sim_in.foff_hz; 103 hf_sim = sim_in.hf_sim; 104 hf_delay = floor(sim_in.hf_delay_ms*Fs/1000); 105 plot_type = sim_in.plot_type; 106 107 % work out gain for HF model 108 % e = sum((g*s)^2) = g*g*sum(s^2) = N, g = sqrt(N/sum(s^2)) 109 % compute so e=N 110 111 s1 = spread(1:frames*M); 112 s2 = [zeros(hf_delay,1); spread_2ms(1:frames*M)]; 113 s2 = s2(1:frames*M); 114 115 p = (s1+s2)'*(s1+s2); 116 hf_gain = sqrt(frames*M/p); 117 p2 = (hf_gain*(s1+s2))'*(hf_gain*(s1+s2)); 118 119 if hf_sim 120 channel_model = "HF"; 121 else 122 channel_model = "AWGN"; 123 end 124 125 % spectrogram states 126 127 spec_states.m = 8*M; 128 spec_states.nfft2 = 2 ^ ceil(log2(spec_states.m/2)); 129 spec_states.dec = 4; 130 spec_states.sample_memory = zeros(1, spec_states.m); 131 spec_states.lower_clip_dB = -30; 132 133 printf("\n%s\n", sim_in.test_name); 134 printf(" Channel EbNo SNR(calc) SNR(meas) SD(Hz) Hits Hits(%%) Result\n"); 135 136 % --------------------------------------------------------------------- 137 % Main loop 138 % --------------------------------------------------------------------- 139 140 for ne = 1:length(EbNovec) 141 EbNo_dB = EbNovec(ne); 142 [Nsd SNR] = calc_Nsd_from_EbNo(EbNo_dB); 143 hits = 0; 144 145 tx_filt = zeros(Nc,M); 146 prev_tx_symbols = ones(Nc+1,1); 147 148 tx_fdm_log = []; 149 rx_fdm_log = []; 150 pilot_lpf1_log = []; 151 S1_log = []; 152 rx_fdm_delay = zeros(M+Ndelay,1); 153 154 % freq offset simulation states 155 156 phase_offset = 1; 157 Nmedian = 20; 158 foff_median=zeros(1,Nmedian); 159 160 % hf sim states 161 162 path2 = zeros(1,hf_delay+M); 163 sum_sig = 0; 164 sum_noise = 0; 165 166 % state machine 167 state = 0; 168 fest_current = 0; 169 fdelta = 5; 170 candidate_thresh = 10; 171 foff_est_thresh_prev = 0; 172 173 for f=1:frames 174 175 % ------------------- Modulator ------------------- 176 177 tx_bits = get_test_bits(Nc*Nb); 178 tx_symbols = bits_to_psk(prev_tx_symbols, tx_bits, 'dqpsk'); 179 180 % simulate BPF filtering of +/- 200 Hz 181 % tx_symbols(1:6) = 0; tx_symbols(9:Nc) = 0; 182 183 prev_tx_symbols = tx_symbols; 184 tx_baseband = tx_filter(tx_symbols); 185 tx_fdm = fdm_upconvert(tx_baseband); 186 tx_fdm_log = [tx_fdm_log real(tx_fdm)]; 187 188 % ------------------- Channel simulation ------------------- 189 190 % frequency offset 191 192 for i=1:M 193 freq_offset = exp(j*2*pi*foff_hz(f)/Fs); 194 phase_offset *= freq_offset; 195 tx_fdm(i) = phase_offset*tx_fdm(i); 196 end 197 198 % optional HF channel sim 199 200 if hf_sim 201 path1 = tx_fdm .* conj(spread(f*M+1:f*M+M)'); 202 203 path2(1:hf_delay) = path2(M+1:hf_delay+M); 204 path2(hf_delay+1:hf_delay+M) = tx_fdm .* conj(spread_2ms(f*M+1:f*M+M)'); 205 206 tx_fdm = hf_gain*(path1 + path2(1:M)); 207 end 208 sum_sig += tx_fdm * tx_fdm'; 209 210 rx_fdm = real(tx_fdm); 211 212 % AWGN noise 213 214 noise = Nsd*randn(1,M); 215 sum_noise += noise * noise'; 216 rx_fdm += noise; 217 rx_fdm_log = [rx_fdm_log rx_fdm]; 218 219 % Fixed Delay 220 221 rx_fdm_delay(1:Ndelay) = rx_fdm_delay(M+1:M+Ndelay); 222 rx_fdm_delay(Ndelay+1:M+Ndelay) = rx_fdm; 223 224 % ------------------- Freq Offset Est ------------------- 225 226 % frequency offset estimation and correction, need to call 227 % rx_est_freq_offset even in track mode to keep states updated 228 229 [pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = ... 230 get_pilot(pilot_lut_index, prev_pilot_lut_index, M); 231 [foff_est S1 S2] = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot, M); 232 pilot_lpf1_log = [pilot_lpf1_log pilot_lpf1(Npilotlpf-M+1:Npilotlpf)]; 233 S1_log(f,:) = fftshift(S1); 234 S2_log(f,:) = fftshift(S2); 235 236 % raw estimate 237 238 foff_log(ne,f) = foff_est; 239 maxS1_log(ne,f) = max(S1.*conj(S1)/(S1*S1')); 240 maxS2_log(ne,f) = max(S2.*conj(S2)/(S2*S2')); 241 242 % median filter post-processed 243 244 foff_median(1:Nmedian-1) = foff_median(2:Nmedian); 245 foff_median(Nmedian) = foff_est; 246 foff_median_log(ne,f) = foff_coarse = median(foff_median); 247 248 % state machine post-processed 249 250 next_state = state; 251 if state == 0 252 if abs(foff_est - fest_current) > fdelta 253 fest_candidate = foff_est; 254 candidate_count = 0; 255 next_state = 1; 256 end 257 end 258 if state == 1 259 if abs(foff_est - fest_candidate) > fdelta 260 next_state = 0; 261 end 262 candidate_count++; 263 if candidate_count > candidate_thresh 264 fest_current = fest_candidate; 265 next_state = 0; 266 end 267 end 268 state = next_state; 269 foff_statemach_log(ne,f) = fest_current; 270 271 % threshold post processed 272 273 if (maxS1_log(ne,f) > 0.06) || (maxS2_log(ne,f) > 0.06) 274 %if (maxS1_log(ne,f) > 0.08) 275 foff_thresh_log(ne,f) = foff_est; 276 else 277 foff_thresh_log(ne,f) = foff_est_thresh_prev; 278 end 279 foff_est_thresh_prev = foff_thresh_log(ne,f); 280 281 % hit/miss stats 282 fest_current = foff_statemach_log(ne,f); 283 if (f > startup_delay) && (abs(fest_current - foff_hz(f)) < allowable_error) 284 hits++; 285 end 286 287 if length(EbNovec) == 1 288 [spectrogram(f,:) spec_states] = spectrogram_update(rx_fdm, M, spec_states); 289 end 290 end 291 292 % results for this EbNo value 293 294 sim_out.foff_sd(ne) = std(foff_log(ne,startup_delay:frames)); 295 sim_out.hits = hits; 296 sim_out.hits_percent = 100*sim_out.hits/(frames-startup_delay); 297 sim_out.SNRvec(ne) = SNR; 298 sim_out.tx_fdm_log = tx_fdm_log; 299 sim_out.rx_fdm_log = rx_fdm_log; 300 301 % noise we have measures is 4000 Hz wide, we want noise in 3000 Hz BW 302 303 snr_meas = 10*log10(sum_sig/(sum_noise*4000/3000)); 304 305 printf(" %6s %5.2f % -5.2f % -5.2f %3.2f %d %3.2f ", ... 306 channel_model, EbNo_dB, SNR, snr_meas, sim_out.foff_sd(ne), sim_out.hits, sim_out.hits_percent); 307 308 if sim_out.hits_percent == 100 309 printf("PASS\n"); 310 else 311 printf("FAIL\n"); 312 figure(5) 313 clf 314 plot(abs(foff_statemach_log(ne,:) - foff_hz < allowable_error)); 315 end 316 317 % plots if single dimension vector 318 319 if length(EbNovec) == 1 320 fmin = -200; fmax = 200; 321 figure(1) 322 clf; 323 subplot(411) 324 plot(foff_log(ne,:)) 325 axis([1 frames fmin fmax]); 326 ylabel("Foff raw") 327 subplot(412) 328 plot(foff_median_log(ne,:)) 329 axis([1 frames fmin fmax]); 330 ylabel("Foff median") 331 subplot(413) 332 plot(foff_statemach_log(ne,:),'g') 333 ylabel("Foff state") 334 axis([1 frames fmin fmax]); 335 subplot(414) 336 plot(foff_thresh_log(ne,:)) 337 ylabel("Foff thresh") 338 axis([1 frames fmin fmax]); 339 xlabel("Frames") 340 grid; 341 342 figure(2) 343 clf; 344 plot(maxS1_log(ne,:)); 345 axis([1 frames 0 0.2]); 346 xlabel("Frames") 347 ylabel("max(abs(S1/S2))") 348 grid; 349 hold on; 350 plot(maxS2_log(ne,:),'g'); 351 hold off; 352 353 figure(3) 354 [n m] = size(S1_log); 355 if strcmp(plot_type,"mesh") 356 mesh(-200+400*(0:m-1)/256,1:n,abs(S1_log(:,:))); 357 xlabel('Freq (Hz)'); ylabel('Frame num'); zlabel("max(abs(S1))") 358 else 359 imagesc(1:n,-200+400*(0:(m-1))/m,abs(S1_log(:,:))'); 360 set(gca,'YDir','normal') 361 ylabel('Freq (Hz)'); xlabel('Frame num'); 362 axis([1 n -200 200]) 363 end 364 365 figure(4) 366 clf 367 [n m] = size(spectrogram); 368 if strcmp(plot_type,"mesh") 369 mesh((4000/m)*(1:m),1:n,spectrogram); 370 xlabel('Freq (Hz)'); ylabel('Frame num'); zlabel('Amplitude (dB)'); 371 else 372 imagesc(1:n,(4000/m)*(1:m),spectrogram') 373 set(gca,'YDir','normal') 374 ylabel('Freq (Hz)'); xlabel('Frame num'); 375 axis([1 n 500 2500]) 376 end 377 378 sim_out.spec = spectrogram; 379 sim_out.tx_fdm_log = spectrogram; 380 end 381 end 382end 383 384% --------------------------------------------------------------------- 385% Run Automated Tests 386% --------------------------------------------------------------------- 387 388more off; 389 390function test1 391 global M; 392 global Rs; 393 394 sim_in.test_name = "Test 1: range of Eb/No (SNRs) in AWGN channel"; 395 sim_in.EbNovec = [3:10 99]; 396 sim_in.delay = M/2; 397 sim_in.frames = Rs*3; 398 sim_in.foff_hz(1:sim_in.frames) = 50; 399 sim_in.startup_delay = 0.5*Rs; 400 sim_in.allowable_error = 5; 401 sim_in.hf_sim = 0; 402 sim_in.hf_delay_ms = 2; 403 sim_in.delay = M/2; 404 sim_in.plot_type = "waterfall"; 405 406 sim_out = freq_off_est_test(sim_in); 407 408 figure(5) 409 clf 410 subplot(211) 411 plot(sim_in.EbNovec, sim_out.foff_sd) 412 hold on; 413 plot(sim_in.EbNovec, sim_out.foff_sd,'+') 414 hold off; 415 xlabel("Eb/No (dB)") 416 ylabel("Std Dev (Hz)") 417 axis([(min(sim_in.EbNovec)-1) (max(sim_in.EbNovec)+1) -1 10]); 418 419 subplot(212) 420 plot(sim_out.SNRvec,sim_out.foff_sd) 421 hold on; 422 plot(sim_out.SNRvec,sim_out.foff_sd,'+') 423 hold off; 424 xlabel("SNR (dB)") 425 ylabel("Std Dev (Hz)") 426 axis([(min(sim_out.SNRvec)-1) (max(sim_out.SNRvec)+1) -1 10]); 427end 428 429 430function test2 431 sim_in.test_name = "Test 2: range of Eb/No (SNRs) in HF multipath channel" 432 sim_in.EbNovec = 0:10; 433 sim_in.delay = 2; 434 sim_in.hf_sim = 1; 435 sim_in.hf_delay_ms = 2; 436 sim_in.frames = Rs*2; 437 sim_in.foff_hz = 0; 438 sim_in.startup_delay = Rs/2; 439 sim_in.allowable_error = 5; 440 441 sim_out = freq_off_est_test(sim_in); 442 443 figure(5) 444 clf 445 subplot(211) 446 plot(sim_in.EbNovec,sim_out.foff_sd) 447 hold on; 448 plot(sim_in.EbNovec,sim_out.foff_sd,'+') 449 hold off; 450 xlabel("Eb/No (dB)") 451 ylabel("Std Dev") 452 axis([(min(sim_in.EbNovec)-1) (max(sim_in.EbNovec)+1) -1 10]); 453end 454 455function test3 456 global M; 457 global Rs; 458 459 sim_in.test_name = "Test 3: 30 Seconds in HF multipath channel at 0dB-ish SNR"; 460 sim_in.EbNovec = 13; 461 sim_in.hf_sim = 0; 462 sim_in.hf_delay_ms = 2; 463 sim_in.delay = M/2; 464 sim_in.frames = Rs; 465 sim_in.foff_hz(1:sim_in.frames) = -50; 466 sim_in.startup_delay = Rs; % allow 1 second in heavily faded channels 467 sim_in.allowable_error = 5; 468 sim_in.plot_type = "mesh"; 469 sim_out = freq_off_est_test(sim_in); 470endfunction 471 472function animated_gif 473 figure(4) 474 for i=5:5:360 475 view(i,45) 476 filename=sprintf('fdmdv_fig%05d.png',i); 477 print(filename); 478 end 479 if 0 480 for i=90:-5:-270 481 view(45,i) 482 filename=sprintf('fdmdv_fig%05d.png',i); 483 print(filename); 484 end 485 end 486endfunction 487 488test3; 489 490
