From 65a278e59a404819df49e30e7e1dcdca87f7a612 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Fri, 13 Apr 2012 06:26:57 +0000 Subject: [PATCH] a little refactoring, getting ready for C port git-svn-id: https://svn.code.sf.net/p/freetel/code@367 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/fdmdv.m | 197 +++++++++++++----------------- codec2-dev/octave/fdmdv_ut.m | 52 +------- codec2-dev/octave/gen_rn_coeffs.m | 40 ++++++ 3 files changed, 129 insertions(+), 160 deletions(-) create mode 100644 codec2-dev/octave/gen_rn_coeffs.m diff --git a/codec2-dev/octave/fdmdv.m b/codec2-dev/octave/fdmdv.m index 7f65353c..5c6c4dab 100644 --- a/codec2-dev/octave/fdmdv.m +++ b/codec2-dev/octave/fdmdv.m @@ -18,49 +18,94 @@ randn('state',1); global Fs = 8000; % sample rate in Hz global T = 1/Fs; % sample period in seconds global Rs = 50; % symbol rate in Hz -global Ts = 1/Rs; % symbol period in seconds global Nc = 14; % number of carriers global Nb = 2; % Bits/symbol for QPSK modulation global Rb = Nc*Rs*Nb; % bit rate global M = Fs/Rs; % oversampling factor global Nsym = 4; % number of symbols to filter over global Fsep = 75; % Separation between carriers (Hz) -global Fcentre = 1200; % Centre frequency, Nc/2 below this, N/c above (Hz) +global Fcentre = 1200; % Centre frequency, Nc/2 carriers below this, N/c carriers above (Hz) global Nt = 5; % number of symbols we estimate timing over global P = 4; % oversample factor used for rx symbol filtering +global Nfilter = Nsym*M; +global Nfiltertiming = M+Nfilter+M; -% Generate root raised cosine (Root Nyquist) filter --------------- +% root raised cosine (Root Nyquist) filter -% thanks http://www.dsplog.com/db-install/wp-content/uploads/2008/05/raised_cosine_filter.m +global gt_alpha5_root = gen_rn_coeffs(alpha, Rs, Nsym, M); -alpha = 0.5; % excess bandwidth -n = -Nsym*Ts/2:T:Nsym*Ts/2; -global Nfilter = Nsym*M; -global Nfiltertiming = M+Nfilter+M; +% Initialise ---------------------------------------------------- -sincNum = sin(pi*n/Ts); % numerator of the sinc function -sincDen = (pi*n/Ts); % denominator of the sinc function -sincDenZero = find(abs(sincDen) < 10^-10); -sincOp = sincNum./sincDen; -sincOp(sincDenZero) = 1; % sin(pix/(pix) =1 for x =0 - -cosNum = cos(alpha*pi*n/Ts); -cosDen = (1-(2*alpha*n/Ts).^2); -cosDenZero = find(abs(cosDen)<10^-10); -cosOp = cosNum./cosDen; -cosOp(cosDenZero) = pi/4; -gt_alpha5 = sincOp.*cosOp; -Nfft = 4096; -GF_alpha5 = fft(gt_alpha5,Nfft)/M; -% sqrt causes stop band to be amplifed, this hack pushes it down again -for i=1:Nfft - if (abs(GF_alpha5(i)) < 0.02) - GF_alpha5(i) *= 0.001; - endif +global pilot_bit; +pilot_bit = 0; % current value of pilot bit + +global tx_filter_memory; +tx_filter_memory = zeros(Nc+1, Nfilter); +global rx_filter_memory; +rx_filter_memory = zeros(Nc+1, Nfilter); + +% phasors used for up and down converters + +global freq; +freq = zeros(Nc+1,1); +for c=1:Nc/2 + carrier_freq = (-Nc/2 - 1 + c)*Fsep + Fcentre; + freq(c) = exp(j*2*pi*carrier_freq/Fs); end -GF_alpha5_root = sqrt(abs(GF_alpha5)) .* exp(j*angle(GF_alpha5)); -ifft_GF_alpha5_root = ifft(GF_alpha5_root); -global gt_alpha5_root = real((ifft_GF_alpha5_root(1:Nfilter))); +for c=Nc/2+1:Nc + carrier_freq = (-Nc/2 + c)*Fsep + Fcentre; + freq(c) = exp(j*2*pi*carrier_freq/Fs); +end + +freq(Nc+1) = exp(j*2*pi*Fcentre/Fs); + +% Spread initial FDM carrier phase out as far as possible. This +% helped PAPR for a few dB. We don't need to adjust rx phase as DQPSK +% takes care of that. + +global phase_tx; +%phase_tx = ones(Nc+1,1); +phase_tx = exp(j*2*pi*(0:Nc)/(Nc+1)); +global phase_rx; +phase_rx = ones(Nc+1,1); + +% Freq offset estimator constants + +global Mpilotfft = 256; +global Npilotcoeff = 30; % number of pilot LPF coeffs +global pilot_coeff = fir1(Npilotcoeff, 200/(Fs/2))'; % 200Hz LPF +global Npilotbaseband = Npilotcoeff + 4*M; % number of pilot baseband samples reqd for pilot LPF +global Npilotlpf = 4*M; % number of samples we DFT pilot over, pilot est window + +% Freq offset estimator states + +global pilot_baseband1; +global pilot_baseband2; +pilot_baseband1 = zeros(1, Npilotbaseband); % pilot baseband samples +pilot_baseband2 = zeros(1, Npilotbaseband); % pilot baseband samples +global pilot_lpf1 +global pilot_lpf2 +pilot_lpf1 = zeros(1, Npilotlpf); % LPF pilot samples +pilot_lpf2 = zeros(1, Npilotlpf); % LPF pilot samples + +% Timing estimator states + +global rx_filter_mem_timing; +rx_filter_mem_timing = zeros(Nc+1, Nt*P); +global rx_baseband_mem_timing; +rx_baseband_mem_timing = zeros(Nc+1, Nfiltertiming); + +% Test bit stream constants + +global Ntest_bits = Nc*Nb*4; % length of test sequence +global test_bits = rand(1,Ntest_bits) > 0.5; + +% Test bit stream state variables + +global current_test_bit = 1; +current_test_bit = 1; +global rx_test_bits_mem; +rx_test_bits_mem = zeros(1,Ntest_bits); % Functions ---------------------------------------------------- @@ -298,7 +343,7 @@ endfunction % LPF and peak pick part of freq est, put in a function as we call it twice -function [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf, nin) +function [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf, nin) global M; global Npilotlpf; global Npilotcoeff; @@ -335,7 +380,7 @@ function [foff imax pilot_lpf] = lpf_peak_pick(pilot_baseband, pilot_lpf, nin) endfunction -% Estimate optimum timing offset, and symbol receive symbols +% Estimate optimum timing offset, re-filter receive symbols function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin) global M; @@ -359,15 +404,12 @@ function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin) % update buffer of Nt rate P filtered symbols rx_filter_mem_timing(:,1:(Nt-1)*P+adjust) = rx_filter_mem_timing(:,P+1-adjust:Nt*P); - %size((Nt-1)*P+1+adjust:Nt*P) - %size(rx_filt) - %adjust - %nin - rx_filter_mem_timing(:,(Nt-1)*P+1+adjust:Nt*P) = rx_filt(1:Nc,:); + rx_filter_mem_timing(:,(Nt-1)*P+1+adjust:Nt*P) = rx_filt(:,:); % sum envelopes of all carriers - env = sum(abs(rx_filter_mem_timing(:,:))); + env = sum(abs(rx_filter_mem_timing(:,:))); % use all Nc+1 carriers for timing + %env = abs(rx_filter_mem_timing(Nc+1,:)); % just use BPSK pilot [n m] = size(env); % The envelope has a frequency component at the symbol rate. The @@ -375,9 +417,9 @@ function [rx_symbols rx_timing] = rx_est_timing(rx_filt, rx_baseband, nin) % single DFT at frequency 2*pi/P x = env * exp(-j*2*pi*(0:m-1)/P)'; - + % map phase to estimated optimum timing instant at rate M - % the M/2 + 1 part was adjusted by experment, I know not why.... + % the M/4 part was adjusted by experiment, I know not why.... rx_timing = angle(x)*M/(2*pi) + M/4; if (rx_timing > M) @@ -519,79 +561,6 @@ function [sync bit_errors] = put_test_bits(rx_bits) endfunction -% Initialise ---------------------------------------------------- - -global pilot_bit; -pilot_bit = 0; % current value of pilot bit - -global tx_filter_memory; -tx_filter_memory = zeros(Nc+1, Nfilter); -global rx_filter_memory; -rx_filter_memory = zeros(Nc+1, Nfilter); - -% phasors used for up and down converters - -global freq; -freq = zeros(Nc+1,1); -for c=1:Nc/2 - carrier_freq = (-Nc/2 - 1 + c)*Fsep + Fcentre; - freq(c) = exp(j*2*pi*carrier_freq/Fs); -end -for c=Nc/2+1:Nc - carrier_freq = (-Nc/2 + c)*Fsep + Fcentre; - freq(c) = exp(j*2*pi*carrier_freq/Fs); -end - -freq(Nc+1) = exp(j*2*pi*Fcentre/Fs); - -% Spread initial FDM carrier phase out as far as possible. This -% helped PAPR for a few dB. We don't need to adjust rx phase as DQPSK -% takes care of that. - -global phase_tx; -%phase_tx = ones(Nc+1,1); -phase_tx = exp(j*2*pi*(0:Nc)/(Nc+1)); -global phase_rx; -phase_rx = ones(Nc+1,1); - -% Freq offset estimator constants - -global Mpilotfft = 256; -global Npilotcoeff = 30; % number of pilot LPF coeffs -global pilot_coeff = fir1(Npilotcoeff, 200/(Fs/2))'; % 200Hz LPF -global Npilotbaseband = Npilotcoeff + 4*M; % number of pilot baseband samples reqd for pilot LPF -global Npilotlpf = 4*M; % number of samples we DFT pilot over, pilot est window - -% Freq offset estimator states constants - -global pilot_baseband1; -global pilot_baseband2; -pilot_baseband1 = zeros(1, Npilotbaseband); % pilot baseband samples -pilot_baseband2 = zeros(1, Npilotbaseband); % pilot baseband samples -global pilot_lpf1 -global pilot_lpf2 -pilot_lpf1 = zeros(1, Npilotlpf); % LPF pilot samples -pilot_lpf2 = zeros(1, Npilotlpf); % LPF pilot samples - -% Timing estimator states - -global rx_filter_mem_timing; -rx_filter_mem_timing = zeros(Nc, Nt*P); -global rx_baseband_mem_timing; -rx_baseband_mem_timing = zeros(Nc+1, Nfiltertiming); - -% Test bit stream constants - -global Ntest_bits = Nc*Nb*4; % length of test sequence -global test_bits = rand(1,Ntest_bits) > 0.5; - -% Test bit stream state variables - -global current_test_bit = 1; -current_test_bit = 1; -global rx_test_bits_mem; -rx_test_bits_mem = zeros(1,Ntest_bits); - % Generate M samples of DBPSK pilot signal for Freq offset estimation @@ -662,7 +631,7 @@ function pilot_lut = generate_pilot_lut pilot_lut = [pilot_lut pilot]; end - % discard first 4 symbols as filer memory is filling, just keep last + % discard first 4 symbols as filter memory is filling, just keep last % four symbols pilot_lut = pilot_lut(4*M+1:M*F); diff --git a/codec2-dev/octave/fdmdv_ut.m b/codec2-dev/octave/fdmdv_ut.m index ad6df497..012f9359 100644 --- a/codec2-dev/octave/fdmdv_ut.m +++ b/codec2-dev/octave/fdmdv_ut.m @@ -13,7 +13,7 @@ fdmdv; % load modem code frames = 100; EbNo_dB = 7.3; -Foff_hz = 10; +Foff_hz = 0; modulation = 'dqpsk'; hpa_clip = 150; @@ -143,9 +143,9 @@ for f=1:frames % 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; + 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 @@ -261,10 +261,10 @@ title('Scatter Diagram'); figure(2) clf; subplot(211) -plot(rx_timing_log(15:m)) +plot(rx_timing_log) title('timing offset (samples)'); subplot(212) -plot(foff_log(15:m)) +plot(foff_log) title('Freq offset (Hz)'); figure(3) @@ -293,43 +293,3 @@ plot(test_frame_sync_log); axis([0 frames 0 1.5]); title('Test Frame Sync') -% TODO -% + handling sample slips, extra plus/minus samples -% + simulating sample clock offsets -% + timing, frequency offset get function -% + memory of recent tx and rx signal for spectrogram -% + scatter diagram get function - -% DPSK -% add phase offset, frequency offset, timing offset -% fading simulator -% file I/O to test with Codec -% code to measure sync recovery time -% dump file type plotting & instrumentation -% determine if error pattern is bursty -% HF channel simulation -% Offset or pi/4 QPSK and tests with real tx HPA -% real time SNR get function -% -% phase estimator not working too well and would need a UW -% to resolve ambiguity. But this is probably worth it for -% 3dB. Test with small freq offset - -% Implementation loss BER issues: -% QPSK mapping -% Single sided noise issues -% interference between carriers due to excess BW -% Crappy RN coeffs -% timing recovery off by one -% Use a true PR sequence -% Sensitivity to Fs -% Try BPSK -% second term of QPSK eqn - -% Faster sync: -% -% 1/ Maybe Ask Bill how we can sync in less than 10 symbols? What to -% put in filter memories? If high SNR should sync very quickly -% Maybe start feeding in symbols to middle of filter memory? Or do timing -% sync before rx filtering at start? - diff --git a/codec2-dev/octave/gen_rn_coeffs.m b/codec2-dev/octave/gen_rn_coeffs.m new file mode 100644 index 00000000..b317d02e --- /dev/null +++ b/codec2-dev/octave/gen_rn_coeffs.m @@ -0,0 +1,40 @@ +% gen_rn_coeffs.m +% David Rowe 13 april 2012 +% +% Generate root raised cosine (Root Nyquist) filter coefficients +% thanks http://www.dsplog.com/db-install/wp-content/uploads/2008/05/raised_cosine_filter.m + +function coeffs = gen_rn_coeffs(alpha, Rs, Nsym, M) + + Ts = 1/Rs; + + n = -Nsym*Ts/2:T:Nsym*Ts/2; + Nfilter = Nsym*M; + Nfiltertiming = M+Nfilter+M; + + sincNum = sin(pi*n/Ts); % numerator of the sinc function + sincDen = (pi*n/Ts); % denominator of the sinc function + sincDenZero = find(abs(sincDen) < 10^-10); + sincOp = sincNum./sincDen; + sincOp(sincDenZero) = 1; % sin(pix/(pix) =1 for x =0 + + cosNum = cos(alpha*pi*n/Ts); + cosDen = (1-(2*alpha*n/Ts).^2); + cosDenZero = find(abs(cosDen)<10^-10); + cosOp = cosNum./cosDen; + cosOp(cosDenZero) = pi/4; + gt_alpha5 = sincOp.*cosOp; + Nfft = 4096; + GF_alpha5 = fft(gt_alpha5,Nfft)/M; + + % sqrt causes stop band to be amplifed, this hack pushes it down again + + for i=1:Nfft + if (abs(GF_alpha5(i)) < 0.02) + GF_alpha5(i) *= 0.001; + endif + end + GF_alpha5_root = sqrt(abs(GF_alpha5)) .* exp(j*angle(GF_alpha5)); + ifft_GF_alpha5_root = ifft(GF_alpha5_root); + coeffs = real((ifft_GF_alpha5_root(1:Nfilter))); +endfunction -- 2.25.1