From 7ab3ed373657679fa6138504501a7e5f1a33b257 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Sat, 22 Dec 2012 06:48:59 +0000 Subject: [PATCH] fist pass at coherent PSK, but not possible I think as PSK symbol symbols for some carriers repeats after just 4 symbols and is not very unique git-svn-id: https://svn.code.sf.net/p/freetel/code@1149 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/fdmdv_ut_coh.m | 412 +++++++++++++++++++++++++++++++ 1 file changed, 412 insertions(+) create mode 100644 codec2-dev/octave/fdmdv_ut_coh.m diff --git a/codec2-dev/octave/fdmdv_ut_coh.m b/codec2-dev/octave/fdmdv_ut_coh.m new file mode 100644 index 00000000..7450b2fd --- /dev/null +++ b/codec2-dev/octave/fdmdv_ut_coh.m @@ -0,0 +1,412 @@ +% fdmdv_ut_coh.m +% + +% Unit Test program for coherent version of FDMDV modem. Used to +% build up the ability to test coherent demodulation of FDMDV +% signals sampled off air. These signals are differentially encoded +% but we can treat the symbols after the diff encoder as PSK symbols. +% +% We keep most of the existing DPSK modem to handle acquisition, frame sync, +% and just the the PSK demo in parallel. The goal here is to measure the BER +% of the test data using coherent PSK, it's not actually a practical modem. + +% Copyright David Rowe 2012 +% This program is distributed under the terms of the GNU General Public License +% Version 2 +% + +fdmdv; % load modem code + +% Simulation Parameters -------------------------------------- + +frames = 100; +EbNo_dB = 73; +Foff_hz = 0; +hpa_clip = 150; + +% ------------------------------------------------------------ + +tx_filt = zeros(Nc,M); +rx_symbols_log = []; +rx_phase_log = 0; +rx_timing_log = 0; +tx_pwr = 0; +noise_pwr = 0; +rx_fdm_log = []; +rx_baseband_log = []; +rx_bits_offset = zeros(Nc*Nb*2); +prev_tx_symbols = ones(Nc+1,1); +prev_rx_symbols = ones(Nc+1,1); +ferr = 0; +foff = 0; +foff_log = []; +tx_baseband_log = []; +tx_fdm_log = []; + +% BER stats + +total_bit_errors = 0; +total_bits = 0; +bit_errors_log = []; +sync_log = []; +test_frame_sync_log = []; +test_frame_sync_state = 0; + +% SNR estimation states + +sig_est = zeros(Nc+1,1); +noise_est = zeros(Nc+1,1); + +% fixed delay simuation + +Ndelay = M+20; +rx_fdm_delay = zeros(Ndelay,1); + +% --------------------------------------------------------------------- +% Eb/No calculations. We need to work out Eb/No for each FDM carrier. +% Total power is sum of power in all FDM carriers +% --------------------------------------------------------------------- + +C = 1; % power of each FDM carrier (energy/sample). Total Carrier power should = Nc*C = Nc +N = 1; % total noise power (energy/sample) of noise source across entire bandwidth + +% Eb = Carrier power * symbol time / (bits/symbol) +% = C *(1/Rs) / 2 +Eb_dB = 10*log10(C) - 10*log10(Rs) - 10*log10(2); + +No_dBHz = Eb_dB - EbNo_dB; + +% Noise power = Noise spectral density * bandwidth +% Noise power = Noise spectral density * Fs/2 for real signals +N_dB = No_dBHz + 10*log10(Fs/2); +Ngain_dB = N_dB - 10*log10(N); +Ngain = 10^(Ngain_dB/20); + +% C/No = Carrier Power/noise spectral density +% = power per carrier*number of carriers / noise spectral density +CNo_dB = 10*log10(C) + 10*log10(Nc) - No_dBHz; + +% SNR in equivalent 3000 Hz SSB channel + +B = 3000; +SNR = CNo_dB - 10*log10(B); + +% freq offset simulation states + +phase_offset = exp(j*0); +freq_offset = exp(j*2*pi*Foff_hz/Fs); +foff_phase = 1; +t = 0; +foff = 0; +fest_state = 0; +track = 0; +track_log = []; + +snr_log = []; + +% PSK demod ---- + +% Set up matrix of tx symbols, one row for each carrier. + +% The differential PSK encoder is driven by a sequence of test bits. +% This test sequence is 4 frames (122 bits long), which means it +% repeats every 4 frames. However the DPSK encoder also has 4 states, +% so it turns out that the transmitted sequence of PSK symbols repeats +% itself every 16 frames. Depending on the test data some carriers +% repeat more often than that, but the sequence as a whole repeats +% every 16 frames. + +% Each row in the matrix is the sequence of symbols for one carrier. +% This repeats itself every 16 symbols. We would like to demodulate +% and measure the BER of this block of 16 symbols. To do this we need +% to correct any phase offset. + +prev_psk_test_symbols = ones(Nc+1,1); +tx_symbols_mem = []; +for i=1:6 + psk_test_symbols = generate_psk_test_symbols(prev_psk_test_symbols); + prev_psk_test_symbols = psk_test_symbols(:,4); + tx_symbols_mem = [tx_symbols_mem psk_test_symbols .* exp(j*pi/4)]; +end + +% PSK BER counters + +psk_total_bits = 0; +psk_total_bit_errors = 0; + +rx_symbols_ph_log = []; + +% --------------------------------------------------------------------- +% Main loop +% --------------------------------------------------------------------- + +for f=1:frames + + % ------------------- + % Modulator + % ------------------- + + tx_bits = get_test_bits(Nc*Nb); + tx_symbols = bits_to_qpsk(prev_tx_symbols, tx_bits, 'dqpsk'); + prev_tx_symbols = tx_symbols; + tx_baseband = tx_filter(tx_symbols); + tx_baseband_log = [tx_baseband_log tx_baseband]; + tx_fdm = fdm_upconvert(tx_baseband); + tx_pwr = 0.9*tx_pwr + 0.1*real(tx_fdm)*real(tx_fdm)'/(M); + + % ------------------- + % Channel simulation + % ------------------- + + % frequency offset + + %Foff_hz += 1/Rs; + Foff = Foff_hz; + for i=1:M + % Time varying freq offset + %Foff = Foff_hz + 100*sin(t*2*pi/(300*Fs)); + %t++; + freq_offset = exp(j*2*pi*Foff/Fs); + phase_offset *= freq_offset; + tx_fdm(i) = phase_offset*tx_fdm(i); + end + + tx_fdm = real(tx_fdm); + + % HPA non-linearity + + tx_fdm(find(abs(tx_fdm) > hpa_clip)) = hpa_clip; + tx_fdm_log = [tx_fdm_log tx_fdm]; + + rx_fdm = tx_fdm; + + % AWGN noise + + noise = Ngain*randn(1,M); + noise_pwr = 0.9*noise_pwr + 0.1*noise*noise'/M; + rx_fdm += noise; + rx_fdm_log = [rx_fdm_log rx_fdm]; + + % Delay + + %rx_fdm_delay(1:Ndelay-M) = rx_fdm_delay(M+1:Ndelay); + %rx_fdm_delay(Ndelay-M+1:Ndelay) = rx_fdm; + rx_fdm_delay = rx_fdm; + + % ------------------- + % Demodulator + % ------------------- + + % frequency offset estimation and correction, need to call + % rx_est_freq_offset even in track mode to keep states updated + + [pilot prev_pilot pilot_lut_index prev_pilot_lut_index] = get_pilot(pilot_lut_index, prev_pilot_lut_index, M); + [foff_course S1 S2] = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot, M); + if track == 0 + foff = foff_course; + end + + foff = 0; % disable for now + + foff_log = [ foff_log foff ]; + foff_rect = exp(j*2*pi*foff/Fs); + + for i=1:M + foff_phase *= foff_rect'; + rx_fdm_delay(i) = rx_fdm_delay(i)*foff_phase; + end + + % baseband processing + + rx_baseband = fdm_downconvert(rx_fdm_delay(1:M), M); + rx_baseband_log = [rx_baseband_log rx_baseband]; + rx_filt = rx_filter(rx_baseband, M); + + [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, M); + rx_timing_log = [rx_timing_log rx_timing]; + + [rx_bits sync foff_fine pd] = qpsk_to_bits(prev_rx_symbols, rx_symbols, 'dqpsk'); + rx_symbols_log = [rx_symbols_log rx_symbols .* exp(j*pi/4)]; + + foff -= 0.5*ferr; + prev_rx_symbols = rx_symbols; + sync_log = [sync_log sync]; + + % freq est state machine + + [track fest_state] = freq_state(sync, fest_state); + track_log = [track_log track]; + + % Update SNR est + + [sig_est noise_est] = snr_update(sig_est, noise_est, pd); + snr_log = [snr_log calc_snr(sig_est, noise_est)]; + + % coherent PSK demod -------------------------------------- + + % update memory of the last Nph rx symbols + + rx_symbols_mem(:,1:Npom-1) = rx_symbols_mem(:,2:Npom); + rx_symbols_mem(:,Npom) = rx_symbols .* exp(j*pi/4); + + % estimate and correct phase of each received symbol in 16 symbol frame + + for c=1:Nc + for s=Nph2+1:Nph2+Nstates % symbols in middle of matrix that we correct + st = s - Nph2; en = s + Nph2; + corr = rx_symbols_mem(c,st:en) * tx_symbols_mem(c,st:en)'; + rx_symbols_mem_ph(c,s) = rx_symbols_mem(c,s) .* exp(-j*angle(corr)); + end + end + + % count bit errors in centre of window + + psk_bit_errors = 0; + psk_bits = Nc*Nstates; + st = Nph2+1; en = Nph2+Nstates; + for c=1:Nc + psk_bit_errors += sum(real(rx_symbols_mem_ph(c,st:en)) .* real(tx_symbols_mem(c,st:en)) < 0); + psk_bit_errors += sum(imag(rx_symbols_mem_ph(c,st:en)) .* imag(tx_symbols_mem(c,st:en)) < 0); + end + printf("psk_bit_errors: %d corr: %f\n", psk_bit_errors, angle(corr)); + ber = psk_bit_errors/psk_bits; + if ber < 0.1 + psk_total_bit_errors += psk_bit_errors; + psk_total_bits += psk_bits; + % phase correction only valid when we are in sync + rx_symbols_ph_log = [rx_symbols_ph_log rx_symbols_mem_ph(:,st:en)]; + end + + % --------------------------------------------------------- + + % count bit errors if we find a test frame + + [test_frame_sync bit_errors] = put_test_bits(test_bits, rx_bits); + + if test_frame_sync == 1 + total_bit_errors = total_bit_errors + bit_errors; + total_bits = total_bits + Ntest_bits; + bit_errors_log = [bit_errors_log bit_errors]; + else + bit_errors_log = [bit_errors_log 0]; + end + + % test frame sync state machine, just for more informative plots + + next_test_frame_sync_state = test_frame_sync_state; + if (test_frame_sync_state == 0) + if (test_frame_sync == 1) + next_test_frame_sync_state = 1; + test_frame_count = 0; + end + end + + if (test_frame_sync_state == 1) + % we only expect another test_frame_sync pulse every 4 symbols + test_frame_count++; + if (test_frame_count == 4) + test_frame_count = 0; + if ((test_frame_sync == 0)) + next_test_frame_sync_state = 0; + end + end + end + test_frame_sync_state = next_test_frame_sync_state; + test_frame_sync_log = [test_frame_sync_log test_frame_sync_state]; +end + +% --------------------------------------------------------------------- +% Print Stats +% --------------------------------------------------------------------- + +ber = total_bit_errors / total_bits; + +% Note Eb/No set point is for Nc data carriers only, excluding pilot. +% This is convenient for testing BER versus Eb/No. Measured Eb/No +% includes power of pilot. Similar for SNR, first number is SNR excluding +% pilot pwr for Eb/No set point, 2nd value is measured SNR which will be a little +% higher as pilot power is included. + +printf("\n"); +printf("Eb/No (meas): %2.2f (%2.2f) dB\n", EbNo_dB, 10*log10(0.25*tx_pwr*Fs/(Rs*Nc*noise_pwr))); +printf("SNR...(meas): %2.2f (%2.2f) dB\n", SNR, calc_snr(sig_est, noise_est)); +printf("\nDPSK\n"); +printf(" bits......: %d\n", total_bits); +printf(" errors....: %d\n", total_bit_errors); +printf(" BER.......: %1.4f\n", ber); +printf("PSK\n"); +printf(" bits......: %d\n", psk_total_bits); +printf(" errors....: %d\n", psk_total_bit_errors); +printf(" BER.......: %1.4f\n", psk_total_bit_errors/psk_total_bits); + +% --------------------------------------------------------------------- +% Plots +% --------------------------------------------------------------------- + +figure(1) +clf; +[n m] = size(rx_symbols_log); +plot(real(rx_symbols_log(1:Nc+1,15:m)),imag(rx_symbols_log(1:Nc+1,15:m)),'+') +%plot(real(rx_symbols_log(2,15:m)),imag(rx_symbols_log(2,15:m)),'+') +axis([-3 3 -3 3]); +title('Scatter Diagram'); + +figure(2) +clf; +subplot(211) +plot(rx_timing_log) +title('timing offset (samples)'); +subplot(212) +plot(foff_log, '-;freq offset;') +hold on; +plot(track_log*75, 'r;course-fine;'); +hold off; +title('Freq offset (Hz)'); + +figure(3) +clf; +subplot(211) +plot(real(tx_fdm_log)); +title('FDM Tx Signal'); +subplot(212) +Nfft=Fs; +S=fft(rx_fdm_log,Nfft); +SdB=20*log10(abs(S)); +plot(SdB(1:Fs/4)) +title('FDM Rx Spectrum'); + +figure(4) +clf; +subplot(311) +stem(sync_log) +axis([0 frames 0 1.5]); +title('BPSK Sync') +subplot(312) +stem(bit_errors_log); +title('Bit Errors for test frames') +subplot(313) +plot(test_frame_sync_log); +axis([0 frames 0 1.5]); +title('Test Frame Sync') + +figure(5) +clf; +subplot(211) +stem(real(rx_symbols_mem(2,:))) +subplot(212) +stem(imag(rx_symbols_mem(2,:))) + +figure(6) +clf; +subplot(211) +stem(real(tx_symbols_mem(2,:))) +subplot(212) +stem(imag(tx_symbols_mem(2,:))) + +figure(7) +clf; +[n m] = size(rx_symbols_ph_log); +plot(real(rx_symbols_ph_log(1:Nc+1,15:m)),imag(rx_symbols_ph_log(1:Nc+1,15:m)),'+') +%plot(real(rx_symbols_ph_log(2,15:m)),imag(rx_symbols_ph_log(2,15:m)),'+') +axis([-3 3 -3 3]); +title('Scatter Diagram - after phase correction'); -- 2.25.1