% LPF cutoff 200Hz, so we can handle max +/- 200 Hz freq offset
pilot_lpf(1:Npilotlpf-nin) = pilot_lpf(nin+1:Npilotlpf);
- j = Npilotcoeff+1;
+ k = Npilotcoeff+1;
for i = Npilotlpf-nin+1:Npilotlpf
- pilot_lpf(i) = pilot_baseband(j-Npilotcoeff+1:j) * pilot_coeff';
- j++;
+ pilot_lpf(i) = pilot_baseband(k-Npilotcoeff+1:k) * pilot_coeff';
+ k++;
end
% decimate to improve DFT resolution, window and DFT
global pilot_baseband2;
global pilot_lpf1;
global pilot_lpf2;
+ global Nbpf;
+ global Nbpfcoeff;
+ global bpf_coeff;
+ global bpf;
+
+if 0
+ % Band Pass filter input so we have (mainly) just the pilots
+
+ bpf(1:Nbpf-nin) = bpf(nin+1:Nbpf);
+ bpf(Nbpf-nin+1:Nbpf) = rx_fdm(1:nin);
+ k = 1;
+ for i = Nbpf-nin+1:Nbpf
+ rx_fdm_bpf(k) = bpf(i-(Nbpfcoeff-1):i) * bpf_coeff';
+ k++;
+ end
+end
% down convert latest nin samples of pilot by multiplying by ideal
% BPSK pilot signal we have generated locally. The peak of the DFT
else
foff = foff2;
end
- foff = foff1;
endfunction
% Freq offset estimator constants
+if 0
+global Nbpfcoeff; % number of input BPF coeffs
+ Nbpfcoeff = 100;
+global bpf_coeff;
+ bpf_coeff = fir1(Nbpfcoeff-1, [Fcentre-100 Fcentre+100]/(Fs/2),'pass')';
+global Nbpf;
+ Nbpf = Nbpfcoeff + M + M/P; % number of pilot baseband samples reqd for pilot LPF
+end
+
global Mpilotfft = 256;
+
global Npilotcoeff; % number of pilot LPF coeffs
- Npilotcoeff = 300; % number of pilot LPF coeffs
+ Npilotcoeff = 30;
global pilot_coeff;
- pilot_coeff = fir1(Npilotcoeff-1, 200/(Fs/2))';% 200Hz LPF
+ pilot_coeff = fir1(Npilotcoeff-1, 100/(Fs/2))';% 200Hz LPF
global Npilotbaseband = Npilotcoeff + M + M/P; % number of pilot baseband samples reqd for pilot LPF
-global Npilotlpf; % number of samples we DFT pilot over, pilot est window
- Npilotlpf = 16*M;
+global Npilotlpf; % number of symbols we DFT pilot over, pilot est window
+ Npilotlpf = 8*M;
% pilot LUT, used for copy of pilot at rx
% Freq offset estimator states
+global bpf;
+ bpf = zeros(1, Nbpf); % BPF pilot input samples
global pilot_baseband1;
global pilot_baseband2;
pilot_baseband1 = zeros(1, Npilotbaseband); % pilot baseband samples
global spread_2ms;
global hf_gain;
- EbNovec = sim_in.EbNovec;
- Ndelay = sim_in.delay;
- frames = sim_in.frames;
+ EbNovec = sim_in.EbNovec;
+ Ndelay = sim_in.delay;
+ frames = sim_in.frames;
startup_delay = sim_in.startup_delay;
allowable_error = sim_in.allowable_error;
- foff_hz = sim_in.foff_hz;
- hf_sim = sim_in.hf_sim;
- hf_delay = floor(sim_in.hf_delay_ms*Fs/1000);
+ foff_hz = sim_in.foff_hz;
+ hf_sim = sim_in.hf_sim;
+ hf_delay = floor(sim_in.hf_delay_ms*Fs/1000);
+ plot_type = sim_in.plot_type;
% work out gain for HF model
% e = sum((g*s)^2) = g*g*sum(s^2) = N, g = sqrt(N/sum(s^2))
hf_gain = sqrt(frames*M/p);
p2 = (hf_gain*(s1+s2))'*(hf_gain*(s1+s2));
+ if hf_sim
+ channel_model = "HF";
+ else
+ channel_model = "AWGN";
+ end
+
% spectrogram states
- spec_states.m = 4*M;
+ spec_states.m = 16*M;
spec_states.nfft2 = 2 ^ ceil(log2(spec_states.m/2));
- spec_states.dec = 4;
+ spec_states.dec = 2;
spec_states.sample_memory = zeros(1, spec_states.m);
spec_states.lower_clip_dB = -30;
+ printf("\n%s\n", sim_in.test_name);
+ printf(" Channel EbNo SNR(calc) SNR(meas) SD(Hz) Hits Hits(%%) Result\n");
+
% ---------------------------------------------------------------------
% Main loop
% ---------------------------------------------------------------------
% freq offset simulation states
phase_offset = 1;
- freq_offset = exp(j*2*pi*foff_hz/Fs);
- foff_phase = 1;
- Nmedian = 25;
+ Nmedian = 20;
foff_median=zeros(1,Nmedian);
% hf sim states
fest_current = 0;
fdelta = 5;
candidate_thresh = 10;
+ foff_est_thresh_prev = 0;
for f=1:frames
tx_bits = get_test_bits(Nc*Nb);
tx_symbols = bits_to_psk(prev_tx_symbols, tx_bits, 'dqpsk');
+
+ % simulate BPF filtering of +/- 200 Hz
+ % tx_symbols(1:6) = 0; tx_symbols(9:Nc) = 0;
+
prev_tx_symbols = tx_symbols;
tx_baseband = tx_filter(tx_symbols);
tx_fdm = fdm_upconvert(tx_baseband);
% frequency offset
- foff = foff_hz;
for i=1:M
- freq_offset = exp(j*2*pi*foff/Fs);
+ freq_offset = exp(j*2*pi*foff_hz(f)/Fs);
phase_offset *= freq_offset;
tx_fdm(i) = phase_offset*tx_fdm(i);
end
[foff_est S1 S2] = rx_est_freq_offset(rx_fdm_delay, pilot, prev_pilot, M);
pilot_lpf1_log = [pilot_lpf1_log pilot_lpf1(Npilotlpf-M+1:Npilotlpf)];
S1_log(f,:) = fftshift(S1);
+ S2_log(f,:) = fftshift(S2);
% raw estimate
foff_log(ne,f) = foff_est;
- maxcf_log(ne,f) = max(abs(S1));
+ maxS1_log(ne,f) = max(S1.*conj(S1)/(S1*S1'));
+ maxS2_log(ne,f) = max(S2.*conj(S2)/(S2*S2'));
% median filter post-processed
if abs(foff_est - fest_candidate) > fdelta
next_state = 0;
end
- candidate_count++;
+ candidate_count++;
if candidate_count > candidate_thresh
fest_current = fest_candidate;
- next_state = 1;
+ next_state = 0;
end
end
state = next_state;
foff_statemach_log(ne,f) = fest_current;
- if (f > startup_delay) && (abs(foff_est - foff_hz) < allowable_error)
+ % threshold post processed
+
+ if (maxS1_log(ne,f) > 0.06) || (maxS2_log(ne,f) > 0.06)
+ %if (maxS1_log(ne,f) > 0.08)
+ foff_thresh_log(ne,f) = foff_est;
+ else
+ foff_thresh_log(ne,f) = foff_est_thresh_prev;
+ end
+ foff_est_thresh_prev = foff_thresh_log(ne,f);
+
+ % hit/miss stats
+ fest_current = foff_statemach_log(ne,f);
+ if (f > startup_delay) && (abs(fest_current - foff_hz(f)) < allowable_error)
hits++;
end
snr_meas = 10*log10(sum_sig/(sum_noise*4000/3000));
- printf("EbNo (dB): %5.2f SNR: % -4.2f % -4.2f std dev (Hz): %3.2f Hits: %d (%3.2f%%)\n", ...
- EbNo_dB, SNR, snr_meas, sim_out.foff_sd(ne), sim_out.hits, sim_out.hits_percent);
+ printf(" %6s %5.2f % -5.2f % -5.2f %3.2f %d %3.2f ", ...
+ channel_model, EbNo_dB, SNR, snr_meas, sim_out.foff_sd(ne), sim_out.hits, sim_out.hits_percent);
+
+ if sim_out.hits_percent == 100
+ printf("PASS\n");
+ else
+ printf("FAIL\n");
+ figure(5)
+ clf
+ plot(abs(foff_statemach_log(ne,:) - foff_hz < allowable_error));
+ end
% plots if single dimension vector
if length(EbNovec) == 1
+ fmin = -200; fmax = 200;
figure(1)
clf;
- subplot(311)
+ subplot(411)
plot(foff_log(ne,:))
- axis([1 frames -200 200]);
- ylabel("Freq offset estimate")
- subplot(312)
+ axis([1 frames fmin fmax]);
+ ylabel("Foff raw")
+ subplot(412)
plot(foff_median_log(ne,:))
- axis([1 frames -200 200]);
- ylabel("Freq offset estimate")
- subplot(313)
- plot(foff_statemach_log(ne,:))
- axis([1 frames -200 200]);
+ axis([1 frames fmin fmax]);
+ ylabel("Foff median")
+ subplot(413)
+ plot(foff_statemach_log(ne,:),'g')
+ ylabel("Foff state")
+ axis([1 frames fmin fmax]);
+ subplot(414)
+ plot(foff_thresh_log(ne,:))
+ ylabel("Foff thresh")
+ axis([1 frames fmin fmax]);
xlabel("Frames")
- ylabel("Freq offset estimate")
grid;
figure(2)
clf;
- plot(maxcf_log(ne,:))
- axis([1 frames -0 300]);
+ plot(maxS1_log(ne,:));
+ axis([1 frames 0 0.2]);
xlabel("Frames")
- ylabel("max(abs(S1))")
+ ylabel("max(abs(S1/S2))")
grid;
+ hold on;
+ plot(maxS2_log(ne,:),'g');
+ hold off;
figure(3)
[n m] = size(S1_log);
- mesh(-200+400*(0:m-1)/256,1:n,abs(S1_log(:,:)));
- xlabel('Freq (Hz)'); ylabel('Frame num'); zlabel("max(abs(S1))")
+ if strcmp(plot_type,"mesh")
+ mesh(-200+400*(0:m-1)/256,1:n,abs(S1_log(:,:)));
+ xlabel('Freq (Hz)'); ylabel('Frame num'); zlabel("max(abs(S1))")
+ else
+ imagesc(1:n,-200+400*(1:(m-1))/m,abs(S1_log(:,:))');
+ set(gca,'YDir','normal')
+ ylabel('Freq (Hz)'); xlabel('Frame num');
+ axis([1 n -200 200])
+ end
figure(4)
clf
[n m] = size(spectrogram);
- lower = floor(500*m/4000); upper = floor(2500*m/4000);
- mesh(lower:upper,1:n,spectrogram(:,lower:upper));
- xlabel('Freq (Hz)'); ylabel('Frame num'); zlabel('Amplitude (dB)');
+ if strcmp(plot_type,"mesh")
+ mesh(1:m,1:n,spectrogram);
+ xlabel('Freq (Hz)'); ylabel('Frame num'); zlabel('Amplitude (dB)');
+ else
+ imagesc(1:n,(4000/m)*(1:m),spectrogram')
+ set(gca,'YDir','normal')
+ ylabel('Freq (Hz)'); xlabel('Frame num');
+ axis([1 n 500 2500])
+ end
sim_out.spec = spectrogram;
sim_out.tx_fdm_log = spectrogram;
% Run Automated Tests
% ---------------------------------------------------------------------
-figure(5); h=freqz(pilot_coeff,1,4000); plot(20*log10(abs(h(1:1000)))); grid
-
more off;
-sim_in.EbNovec = 3;
-sim_in.hf_sim = 1;
-sim_in.hf_delay_ms = 2;
-sim_in.delay = M/2;
-sim_in.frames = Rs*10;
-sim_in.foff_hz = 50;
-sim_in.startup_delay = 10;
-sim_in.allowable_error = 5;
-
-sim_out = freq_off_est_test(sim_in);
-
-% Test 1 - range of Eb/No (SNRs) in AWGN channel
-
function test1
- sim_in.EbNovec = 0:10;
+ global M;
+ global Rs;
+
+ sim_in.test_name = "Test 1: range of Eb/No (SNRs) in AWGN channel";
+ sim_in.EbNovec = [3:10 99];
sim_in.delay = M/2;
- sim_in.frames = 20;
- sim_in.foff_hz = 50;
- sim_in.startup_delay = 10;
+ sim_in.frames = Rs*3;
+ sim_in.foff_hz(1:sim_in.frames) = 50;
+ sim_in.startup_delay = 0.5*Rs;
sim_in.allowable_error = 5;
+ sim_in.hf_sim = 0;
+ sim_in.hf_delay_ms = 2;
+ sim_in.delay = M/2;
+ sim_in.plot_type = "waterfall";
sim_out = freq_off_est_test(sim_in);
- figure(4)
+ figure(5)
clf
subplot(211)
- plot(sim_in.EbNovec,sim_out.foff_sd)
+ plot(sim_in.EbNovec, sim_out.foff_sd)
hold on;
- plot(sim_in.EbNovec,sim_out.foff_sd,'+')
+ plot(sim_in.EbNovec, sim_out.foff_sd,'+')
hold off;
xlabel("Eb/No (dB)")
- ylabel("Std Dev")
+ ylabel("Std Dev (Hz)")
axis([(min(sim_in.EbNovec)-1) (max(sim_in.EbNovec)+1) -1 10]);
subplot(212)
plot(sim_out.SNRvec,sim_out.foff_sd,'+')
hold off;
xlabel("SNR (dB)")
- ylabel("Std Dev")
+ ylabel("Std Dev (Hz)")
axis([(min(sim_out.SNRvec)-1) (max(sim_out.SNRvec)+1) -1 10]);
end
-% Test 2 - range of Eb/No (SNRs) in multipath channel
function test2
+ sim_in.test_name = "Test 2: range of Eb/No (SNRs) in HF multipath channel"
sim_in.EbNovec = 0:10;
sim_in.delay = 2;
- sim_in.hf_sim = 0;
+ sim_in.hf_sim = 1;
sim_in.hf_delay_ms = 2;
- sim_in.frames = Rs*10;
+ sim_in.frames = Rs*2;
sim_in.foff_hz = 0;
- sim_in.startup_delay = 10;
+ sim_in.startup_delay = Rs/2;
sim_in.allowable_error = 5;
sim_out = freq_off_est_test(sim_in);
axis([(min(sim_in.EbNovec)-1) (max(sim_in.EbNovec)+1) -1 10]);
end
+function test3
+ global M;
+ global Rs;
+
+ sim_in.test_name = "Test 3: 30 Seconds in HF multipath channel at 0dB-ish SNR";
+ sim_in.EbNovec = 3;
+ sim_in.hf_sim = 0;
+ sim_in.hf_delay_ms = 2;
+ sim_in.delay = M/2;
+ sim_in.frames = Rs*5;
+ sim_in.foff_hz(1:sim_in.frames) = -50;
+ sim_in.startup_delay = Rs; % allow 1 second in heavily faded channels
+ sim_in.allowable_error = 5;
+ sim_in.plot_type = "waterfall";
+ sim_out = freq_off_est_test(sim_in);
+endfunction
+
+test3;
+
projects / freetel-svn-tracking.git / commitdiff