1% fdmdv_ut.m 2% 3% Unit Test program for FDMDV modem. Useful for general development as it has 4% both tx and rx sides, and basic AWGN channel simulation. 5% 6% Copyright David Rowe 2012 7% This program is distributed under the terms of the GNU General Public License 8% Version 2 9% 10 11fdmdv; % load modem code 12 13% Simulation Parameters -------------------------------------- 14 15frames = 100; 16EbNo_dB = 6.3; 17Foff_hz = -100; 18modulation = 'dqpsk'; 19hpa_clip = 150; 20 21% ------------------------------------------------------------ 22 23more off; 24tx_filt = zeros(Nc,M); 25rx_symbols_log = []; 26rx_phase_log = 0; 27rx_timing_log = 0; 28tx_pwr = 0; 29noise_pwr = 0; 30rx_fdm_log = []; 31rx_baseband_log = []; 32rx_bits_offset = zeros(Nc*Nb*2); 33prev_tx_symbols = ones(Nc+1,1); prev_tx_symbols(Nc+1) = 2; 34prev_rx_symbols = ones(Nc+1,1); 35ferr = 0; 36foff = 0; 37foff_log = []; 38tx_baseband_log = []; 39tx_fdm_log = []; 40 41% BER stats 42 43total_bit_errors = 0; 44total_bits = 0; 45bit_errors_log = []; 46sync_bit_log = []; 47test_frame_sync_log = []; 48test_frame_sync_state = 0; 49 50% SNR estimation states 51 52sig_est = zeros(Nc+1,1); 53noise_est = zeros(Nc+1,1); 54 55% fixed delay simuation 56 57Ndelay = M+20; 58rx_fdm_delay = zeros(Ndelay,1); 59 60% --------------------------------------------------------------------- 61% Eb/No calculations. We need to work out Eb/No for each FDM carrier. 62% Total power is sum of power in all FDM carriers 63% --------------------------------------------------------------------- 64 65C = 1; % power of each FDM carrier (energy/sample). Total Carrier power should = Nc*C = Nc 66N = 1; % total noise power (energy/sample) of noise source across entire bandwidth 67 68% Eb = Carrier power * symbol time / (bits/symbol) 69% = C *(1/Rs) / Nb 70Eb_dB = 10*log10(C) - 10*log10(Rs) - 10*log10(Nb); 71 72No_dBHz = Eb_dB - EbNo_dB; 73 74% Noise power = Noise spectral density * bandwidth 75% Noise power = Noise spectral density * Fs/2 for real signals 76N_dB = No_dBHz + 10*log10(Fs/2); 77Ngain_dB = N_dB - 10*log10(N); 78Ngain = 10^(Ngain_dB/20); 79 80% C/No = Carrier Power/noise spectral density 81% = power per carrier*number of carriers / noise spectral density 82CNo_dB = 10*log10(C) + 10*log10(Nc) - No_dBHz; 83 84% SNR in equivalent 3000 Hz SSB channel 85 86B = 3000; 87SNR = CNo_dB - 10*log10(B); 88 89% freq offset simulation states 90 91phase_offset = 1; 92freq_offset = exp(j*2*pi*Foff_hz/Fs); 93foff_phase = 1; 94t = 0; 95foff = 0; 96fest_state = 0; 97fest_timer = 0; 98sync_mem = zeros(1,Nsync_mem); 99sync = 0; 100sync_log = []; 101 102snr_log = []; 103 104Nspec=1024; 105spec_mem=zeros(1,Nspec); 106SdB = zeros(1,Nspec); 107 108% --------------------------------------------------------------------- 109% Main loop 110% --------------------------------------------------------------------- 111 112for f=1:frames 113 114 % ------------------- 115 % Modulator 116 % ------------------- 117 118 tx_bits = get_test_bits(Nc*Nb); 119 tx_symbols = bits_to_psk(prev_tx_symbols, tx_bits, modulation); 120 prev_tx_symbols = tx_symbols; 121 tx_baseband = tx_filter(tx_symbols); 122 tx_baseband_log = [tx_baseband_log tx_baseband]; 123 tx_fdm = fdm_upconvert(tx_baseband); 124 tx_pwr = 0.9*tx_pwr + 0.1*real(tx_fdm)*real(tx_fdm)'/(M); 125 126 % ------------------- 127 % Channel simulation 128 % ------------------- 129 130 % frequency offset 131 132 %Foff_hz += 1/Rs; 133 Foff = Foff_hz; 134 for i=1:M 135 % Time varying freq offset 136 %Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs)); 137 %t++; 138 freq_offset = exp(j*2*pi*Foff/Fs); 139 phase_offset *= freq_offset; 140 tx_fdm(i) = phase_offset*tx_fdm(i); 141 end 142 143 tx_fdm = real(tx_fdm); 144 145 % HPA non-linearity 146 147 tx_fdm(find(abs(tx_fdm) > hpa_clip)) = hpa_clip; 148 tx_fdm_log = [tx_fdm_log tx_fdm]; 149 150 rx_fdm = tx_fdm; 151 152 % AWGN noise 153 154 noise = Ngain*randn(1,M); 155 noise_pwr = 0.9*noise_pwr + 0.1*noise*noise'/M; 156 rx_fdm += noise; 157 rx_fdm_log = [rx_fdm_log rx_fdm]; 158 159 % update spectrum 160 161 l=length(rx_fdm); 162 spec_mem(1:Nspec-l) = spec_mem(l+1:Nspec); 163 spec_mem(Nspec-l+1:Nspec) = rx_fdm; 164 S=fft(spec_mem.*hanning(Nspec)',Nspec); 165 SdB = 0.9*SdB + 0.1*20*log10(abs(S)); 166 167 168 % ------------------- 169 % Demodulator 170 % ------------------- 171 172 % shift down to complex baseband 173 174 for i=1:M 175 fbb_phase_rx = fbb_phase_rx*fbb_rect'; 176 rx_fdm(i) = rx_fdm(i)*fbb_phase_rx; 177 end 178 mag = abs(fbb_phase_rx); 179 fbb_phase_rx /= mag; 180 181 % frequency offset estimation and correction, need to call rx_est_freq_offset even in sync 182 % mode to keep states updated 183 184 [pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pilot_lut_index, prev_pilot_lut_index, M); 185 [foff_coarse S1 S2] = rx_est_freq_offset(rx_fdm, pilot, prev_pilot, M, !sync); 186 187 if sync == 0 188 foff = foff_coarse; 189 end 190 191 foff_log = [ foff_log foff ]; 192 foff_rect = exp(j*2*pi*foff/Fs); 193 194 for i=1:M 195 foff_phase *= foff_rect'; 196 rx_fdm(i) = rx_fdm(i)*foff_phase; 197 end 198 199 rx_fdm_filter = rxdec_filter(rx_fdm, M); 200 rx_filt = down_convert_and_rx_filter(rx_fdm_filter, M, M/Q); 201 202 [rx_symbols rx_timing] = rx_est_timing(rx_filt, M); 203 rx_timing_log = [rx_timing_log rx_timing]; 204 205 %rx_phase = rx_est_phase(rx_symbols); 206 %rx_phase_log = [rx_phase_log rx_phase]; 207 %rx_symbols = rx_symbols*exp(j*rx_phase); 208 209 [rx_bits sync_bit foff_fine pd] = psk_to_bits(prev_rx_symbols, rx_symbols, modulation); 210 if strcmp(modulation,'dqpsk') 211 rx_symbols_log = [rx_symbols_log pd]; 212 else 213 rx_symbols_log = [rx_symbols_log rx_symbols]; 214 endif 215 foff -= 0.5*foff_fine; 216 217 prev_rx_symbols = rx_symbols; 218 sync_bit_log = [sync_bit_log sync_bit]; 219 220 % freq est state machine 221 222 [sync reliable_sync_bit fest_state fest_timer sync_mem] = freq_state(sync_bit, fest_state, fest_timer, sync_mem); 223 sync_log = [sync_log sync]; 224 225 % Update SNR est 226 227 [sig_est noise_est] = snr_update(sig_est, noise_est, pd); 228 snr_log = [snr_log calc_snr(sig_est, noise_est)]; 229 230 % count bit errors if we find a test frame 231 % Allow 15 frames for filter memories to fill and time est to settle 232 233 [test_frame_sync bit_errors] = put_test_bits(test_bits, rx_bits); 234 235 if test_frame_sync == 1 236 total_bit_errors = total_bit_errors + bit_errors; 237 total_bits = total_bits + Ntest_bits; 238 bit_errors_log = [bit_errors_log bit_errors]; 239 else 240 bit_errors_log = [bit_errors_log 0]; 241 end 242 243 % test frame sync state machine, just for more informative plots 244 245 next_test_frame_sync_state = test_frame_sync_state; 246 if (test_frame_sync_state == 0) 247 if (test_frame_sync == 1) 248 next_test_frame_sync_state = 1; 249 test_frame_count = 0; 250 end 251 end 252 253 if (test_frame_sync_state == 1) 254 % we only expect another test_frame_sync pulse every 4 symbols 255 test_frame_count++; 256 if (test_frame_count == 4) 257 test_frame_count = 0; 258 if ((test_frame_sync == 0)) 259 next_test_frame_sync_state = 0; 260 end 261 end 262 end 263 test_frame_sync_state = next_test_frame_sync_state; 264 test_frame_sync_log = [test_frame_sync_log test_frame_sync_state]; 265end 266 267% --------------------------------------------------------------------- 268% Print Stats 269% --------------------------------------------------------------------- 270 271ber = total_bit_errors / total_bits; 272 273% Peak to Average Power Ratio calcs from http://www.dsplog.com 274 275papr = max(tx_fdm_log.*conj(tx_fdm_log)) / mean(tx_fdm_log.*conj(tx_fdm_log)); 276papr_dB = 10*log10(papr); 277 278% Note Eb/No set point is for Nc data carriers only, excluding pilot. 279% This is convenient for testing BER versus Eb/No. Measured SNR & 280% Eb/No includes power of pilot. Similar for SNR, first number is SNR 281% excluding pilot pwr for Eb/No set point, 2nd value is measured SNR 282% which will be a little higher as pilot power is included. Note current SNR 283% est algorithm only works for QPSK, gives silly values for 8PSK. 284 285printf("Bits/symbol.: %d\n", Nb); 286printf("Num carriers: %d\n", Nc); 287printf("Bit Rate....: %d bits/s\n", Rb); 288printf("Eb/No (meas): %2.2f (%2.2f) dB\n", EbNo_dB, 10*log10(0.25*tx_pwr*Fs/(Rs*Nc*noise_pwr))); 289printf("bits........: %d\n", total_bits); 290printf("errors......: %d\n", total_bit_errors); 291printf("BER.........: %1.4f\n", ber); 292printf("PAPR........: %1.2f dB\n", papr_dB); 293printf("SNR...(meas): %2.2f (%2.2f) dB\n", SNR, calc_snr(sig_est, noise_est)); 294 295% --------------------------------------------------------------------- 296% Plots 297% --------------------------------------------------------------------- 298 299figure(1) 300clf; 301[n m] = size(rx_symbols_log); 302plot(real(rx_symbols_log(1:Nc+1,15:m)),imag(rx_symbols_log(1:Nc+1,15:m)),'+') 303axis([-3 3 -3 3]); 304title('Scatter Diagram'); 305 306figure(2) 307clf; 308subplot(211) 309plot(rx_timing_log) 310title('timing offset'); 311subplot(212) 312plot(foff_log, '-;freq offset;') 313hold on; 314plot(sync_log*75, 'r;Sync State & course(0) fine(1) freq tracking;'); 315hold off; 316title('Freq offset (Hz)'); 317 318figure(3) 319clf; 320subplot(211) 321plot(real(tx_fdm_log)); 322title('FDM Tx Signal'); 323subplot(212) 324plot((0:Nspec/2-1)*Fs/Nspec, SdB(1:Nspec/2) - 20*log10(Nspec/2)) 325axis([0 Fs/2 -40 0]) 326grid 327title('FDM Rx Spectrum'); 328 329figure(4) 330clf; 331subplot(311) 332stem(sync_bit_log) 333axis([0 frames 0 1.5]); 334title('BPSK Sync') 335subplot(312) 336stem(bit_errors_log); 337title('Bit Errors for test frames') 338subplot(313) 339plot(test_frame_sync_log); 340axis([0 frames 0 1.5]); 341title('Test Frame Sync') 342 343figure(5) 344clf 345subplot(211) 346plot(snr_log) 347subplot(212) 348%plot(20*log10(sig_est(1:Nc))-20*log10(sig_est(Nc+1))+6) 349%axis([1 Nc -6 6]); 350sdB_pc = 20*log10(sig_est(1:Nc+1)); 351bar(sdB_pc(1:Nc) - mean(sdB_pc(1:Nc))) 352axis([0 Nc+1 -3 3]); 353