From 1b55d93bc58d8752e3854cb14af7d31a30a85a57 Mon Sep 17 00:00:00 2001 From: drowe67 Date: Thu, 24 Mar 2016 10:16:34 +0000 Subject: [PATCH] cleaned up code, nice compas plot git-svn-id: https://svn.code.sf.net/p/freetel/code@2759 01035d8c-6547-0410-b346-abe4f91aad63 --- codec2-dev/octave/df_mixer.m | 278 +++++++++++++++++++---------------- 1 file changed, 152 insertions(+), 126 deletions(-) diff --git a/codec2-dev/octave/df_mixer.m b/codec2-dev/octave/df_mixer.m index f103573b..51e463c2 100644 --- a/codec2-dev/octave/df_mixer.m +++ b/codec2-dev/octave/df_mixer.m @@ -6,17 +6,15 @@ % antennas. This simulation tests estimation of phase and freq of % mixer LO. % -% [X] any block size, dft size, Fs -% [X] stats for repeated tests -% [X] add modulation -% [ ] sample by sample approach -% + shooting for some processing gain -% [ ] tuning accuracy -% [ ] imperfect mixer with carrier feedthru -% [ ] try with real signal +% [ ] built in hackrf_dc +% [ ] test mode +% [ ] compass needle +% [ ] real time operation (or close to it) fm; +% Generate a FM modulated signal for testing ---------------------------------- + function tx = fm_mod(Fs, fmod) fm_states.Fs = Fs; fm_states.fm_max = 3E3; @@ -36,12 +34,60 @@ function tx = fm_mod(Fs, fmod) endfunction +% Check Filters ------------------------------------------------------------- + +% Run some test signals through all three filters. Real signals used +% for convenience, also works with complex. If another set of filter +% is used you need to prove net reponse is a delay. Fig 2 +% subplot(211) must be a straight line. + +function check_filters(bcentre, Fs, f, bw, w, wlo, t) + x = cos(w*t) + cos((w-wlo)*t) + cos((w-wlo)*t); + ycenter = filter(bcentre,1,x); + ylow = exp(-j*wlo*t) .* filter(bcentre, 1, x .* exp(j*wlo*t)); + yhigh = exp(j*wlo*t) .* filter(bcentre, 1, x .* exp(-j*wlo*t)); + y = ycenter + ylow + yhigh; + + figure(1) + clf + h=freqz(bcentre,1,Fs/2); + plot(20*log10(abs(h))) + hold on; + plot([f-bw/2 f-bw/2 f+bw/2 f+bw/2],[-60 -3 -3 -60],'r'); + hold off; + axis([f-bw f+bw -60 0]) + title('Centre BPF'); + grid; + + figure(2) + clf; + Ns = 200; + Np = 200; + + % y should be a delayed version of x + + subplot(211) + plot(x(Ns:Ns+Np-1)) + hold on; + plot(y(Ns:Ns+Np-1),'g') + hold off; + title('Filter Check'); + + % should see a straight line on this plot + + subplot(212) + plot(x(Ns:Ns+Np-1), y(Ns:Ns+Np-1)) +end + + % Constants ------------------------------------------------------------------- randn('state',1); more off; format short +simulate = 0; % set to 0 for HackRF input, 1 for simulated input + f = 48E3; % signal frequency flo = 32.768E3; % LO frequency m = 10^(-0/20); % modulation index of our mixer, how far each @@ -56,6 +102,7 @@ theta = 0; % phase term constant to ant1 and ant2 alpha = 0; % LO phase Fs = 200E3; % sample rate +Fshrf = 10E6; % HackRF sample rate CNodB = 100; % C/No of carrier, which dominates power % Around C/No = 50dB is min for FM (12dB-ish SINAD) fmod = 0; % Audio modulation on FM signal (0 for carrier) @@ -68,15 +115,10 @@ N = Fs/2; % processing block size CNo = 10^(CNodB/10); var = Fs/CNo; -t=0:N-1; - w = 2*pi*f/Fs; wlo = 2*pi*flo/Fs; wbw = 2*pi*bw/Fs; - -% Check Filters ------------------------------------------------------------- - % FIR filter BPF for centre signal. Use a 15kHz cut off to catch % energy from FM signals. Important: we need identical phase response % for each BPF. I couldn't work out how to design these filter using @@ -85,75 +127,43 @@ wbw = 2*pi*bw/Fs; bcentre = fir2(200, [0, f-bw/2, f, f+bw/2, Fs/2]/(Fs/2), [0 0 1 0 0]); -% Run some test signals through all three filters. Real signals used -% for convenience, also works with complex. If another set of filter -% is used you need to prove net reponse is a delay. Fig 2 -% subplot(211) must be a straight line. - -x = cos(w*t) + cos((w-wlo)*t) + cos((w-wlo)*t); -ycenter = filter(bcentre,1,x); -ylow = exp(-j*wlo*t) .* filter(bcentre, 1, x .* exp(j*wlo*t)); -yhigh = exp(j*wlo*t) .* filter(bcentre, 1, x .* exp(-j*wlo*t)); -y = ycenter + ylow + yhigh; - - -figure(1) -clf -h=freqz(bcentre,1,Fs/2); -plot(20*log10(abs(h))) -hold on; -plot([f-bw/2 f-bw/2 f+bw/2 f+bw/2],[-60 -3 -3 -60],'r'); -hold off; -axis([f-bw f+bw -60 0]) -title('Centre BPF'); -grid; +%check_filters(bcentre, Fs, f, bw, w, wlo, t); xx -figure(2) -clf; -Ns = 200; -Np = 200; - -% y should be a delayed version of x - -subplot(211) -plot(x(Ns:Ns+Np-1)) -hold on; -plot(y(Ns:Ns+Np-1),'g') -hold off; -title('Filter Check'); +% Main ------------------------------------------------------------------------- -% should see a straight line on this plot +if simulate + % Simulate rx signal --------------------------------------------------------- -subplot(212) -plot(x(Ns:Ns+Np-1), y(Ns:Ns+Np-1)) + t=0:N-1; + % tx signal -% Main ------------------------------------------------------------------------- + tx = fm_mod(Fs, fmod); -% tx signal + % simulate rx signals at antenna 1 and 2 by adding noise. Note signal + % at antenna 2 is phase shifted phi due to path difference. Both signals + % experience a time of flight path different theta. -tx = fm_mod(Fs, fmod); + n = sqrt(var/2)*randn(1,N) + j*sqrt(var/2)*randn(1,N); + ant1 = tx*exp(j*(theta+phi)) + n; + ant2 = tx*exp(j*(theta)) + n; -% simulate rx signals at antenna 1 and 2 by adding noise. Note signal -% at antenna 2 is phase shifted phi due to path difference. Both signals -% experience a time of flight path different theta. + % ant2 passed through mixer with local osc frequency wlo and phase alpha -n = sqrt(var/2)*randn(1,N) + j*sqrt(var/2)*randn(1,N); -ant1 = tx*exp(j*(theta+phi)) + n; -ant2 = tx*exp(j*(theta)) + n; + ant2_mix = 2.0*m*(ant2 .* cos(t*wlo + alpha)); -% ant2 passed through mixer with local osc frequency wlo and phase alpha + % The SDR receives the sum of the two signals -ant2_mix = 2.0*m*(ant2 .* cos(t*wlo + alpha)); + rx = ant1 + ant2_mix; -% The SDR receives the sum of the two signals +else + % Off air signal from HackRF ------------------------------------------------------- -rx = ant1 + ant2_mix; + s1 = load_hackrf("df1.iq"); + rx = downsample(rot90(s1), Fshrf/Fs)/127; + t = 1:length(rx); +end -% Uncomment to use HackRF file as input, ow test sig generated baove is used - -rx=hackrf_dc("df1.iq"); rx = rx(1000:length(rx)); % discard initial transients -t = 1:length(rx); Rx = (1/N)*fft(rx,N); % BPF each signal @@ -162,71 +172,87 @@ sam1 = filter(bcentre,1,rx); sam2 = exp( j*wlo*t) .* filter(bcentre, 1, rx .* exp(-j*wlo*t)); sam3 = exp(-j*wlo*t) .* filter(bcentre, 1, rx .* exp( j*wlo*t)); -%sam1 = sign(sam1); sam2 = sign(sam2); sam3 = sign(sam3); - % Do the math on phases using rect math two_phi_est_rect = conj(sam2.*sam3) .* (sam1.^2); phi_est = angle(two_phi_est_rect)/2; phi_est_av = angle(mean(two_phi_est_rect))/2; + printf("phi_est : %3.2f rads %3.1f degrees\n", phi_est_av, phi_est_av*180/pi); % some plots .... -figure(2); -clf; -plot((1:N)*Fs/N, 20*log10(abs(Rx)),'markersize', 10, 'linewidth', 2); -axis([1 Fs -80 0]) -title('Rx signal at SDR input'); - -figure(3) -clf -Sam1 = (1/N)*fft(sam1,N); -Sam2 = (1/N)*fft(sam2,N); -Sam3 = (1/N)*fft(sam3,N); -subplot(311) -plot((1:N)*Fs/N, 20*log10(abs(Sam1)),'markersize', 10, 'linewidth', 2); -axis([1 Fs/2 -60 0]) -title('Band Pass filtered signals'); -subplot(312) -plot((1:N)*Fs/N, 20*log10(abs(Sam2)),'markersize', 10, 'linewidth', 2); -axis([1 Fs/2 -60 0]) -subplot(313) -plot((1:N)*Fs/N, 20*log10(abs(Sam3)),'markersize', 10, 'linewidth', 2); -axis([1 Fs/2 -60 0]) - -figure(4) -Np = 1000; -subplot(311) -plot(real(sam1(1:Np))) -title('Band Pass filtered signals'); -subplot(312) -plot(real(sam2(1:Np))) -subplot(313) -plot(real(sam3(1:Np))) - -figure(5) -clf; -plot(two_phi_est_rect) -xmax = 2*mean(abs(two_phi_est_rect)); -axis([-xmax xmax -xmax xmax]) -title('Scatter Plot'); - -figure(6) -clf; -subplot(211) -plot(real(two_phi_est_rect(1:Np))) -subplot(212) -plot(imag(two_phi_est_rect(1:Np))) -title('two phi est rect') - -figure(7) -clf -hist([angle(two_phi_est_rect)/2 -pi/2 pi/2],100) - -figure(8) -clf -plot(phi_est) -axis([0 length(phi_est(1:Np)) -pi/2 pi/2]); -title('phi est') +figs = 0x1 + 0x20; + +if bitand(figs, 0x1) + figure(1); + clf; + plot((1:N)*Fs/N, 20*log10(abs(Rx)),'markersize', 10, 'linewidth', 2); + axis([1 Fs -80 0]) + title('Rx signal at SDR input'); +end + +if bitand(figs, 0x2) + figure(2) + clf + Sam1 = (1/N)*fft(sam1,N); + Sam2 = (1/N)*fft(sam2,N); + Sam3 = (1/N)*fft(sam3,N); + subplot(311) + plot((1:N)*Fs/N, 20*log10(abs(Sam1)),'markersize', 10, 'linewidth', 2); + axis([1 Fs/2 -60 0]) + title('Band Pass filtered signals'); + subplot(312) + plot((1:N)*Fs/N, 20*log10(abs(Sam2)),'markersize', 10, 'linewidth', 2); + axis([1 Fs/2 -60 0]) + subplot(313) + plot((1:N)*Fs/N, 20*log10(abs(Sam3)),'markersize', 10, 'linewidth', 2); + axis([1 Fs/2 -60 0]) +end + +if bitand(figs, 0x4) + figure(3) + Np = 1000; + subplot(311) + plot(real(sam1(1:Np))) + title('Band Pass filtered signals'); + subplot(312) + plot(real(sam2(1:Np))) + subplot(313) + plot(real(sam3(1:Np))) +end + +if bitand(figs, 0x8) + figure(4) + clf; + hist2d([real(two_phi_est_rect)' imag(two_phi_est_rect)'],50) + xmax = 2*mean(abs(two_phi_est_rect)); + axis([-xmax xmax -xmax xmax]) + title('Scatter Plot'); +end + +if bitand(figs, 0x10) + figure(5) + clf; + subplot(211) + plot(real(two_phi_est_rect(1:Np))) + subplot(212) + plot(imag(two_phi_est_rect(1:Np))) + title('two phi est rect') +end + +if bitand(figs, 0x20) + figure(6) + clf + [r theta] = hist([angle(two_phi_est_rect)/2 -pi/2 pi/2],100); + polar([theta pi+theta],[r r]) +end + +if bitand(figs, 0x40) + figure(7) + clf + plot(phi_est) + axis([0 length(phi_est(1:Np)) -pi/2 pi/2]); + title('phi est') +end -- 2.25.1