doppler spreading functions for HF channel simulation

git-svn-id: https://svn.code.sf.net/p/freetel/code@2652 01035d8c-6547-0410-b346-abe4f91aad63

diff --git a/codec2-dev/octave/doppler_spread.m b/codec2-dev/octave/doppler_spread.m
new file mode 100644 (file)
index 0000000..4bc6c62
--- /dev/null
+++ b/codec2-dev/octave/doppler_spread.m
@@ -0,0 +1,37 @@
+% doppler_spread.m
+% David Rowe Jan 2016
+%
+% Returns gausssian filtered doppler spreading function samples for HF channel
+% modelling.  Used PathSim technical guide as a reference - thanks Moe!
+
+function [spread_FsHz states] = doppler_spread(dopplerSpreadHz, FsHz, Nsam)
+
+  % start with low Fs so we can work with a reasonable filter length
+
+  sigma = dopplerSpreadHz/2;
+  lowFs = 10*dopplerSpreadHz;
+  Ntaps = 100;
+  Nsam_low = Nsam*lowFs/FsHz + Ntaps; % fill filter memory
+
+  % generate gaussian freq response and design filter
+
+  x = 0:0.1:lowFs/2;
+  y = (1/(sigma*sqrt(2*pi)))*exp(-(x.^2)/(2*sigma*sigma));
+  b = fir2(Ntaps-1, x/(lowFs/2), y);
+
+  % generate the spreading samples
+
+  spread_lowFs = filter(b,1,randn(1,Nsam_low) + j*randn(1,Nsam_low));
+
+  % resample to FsHz, scaling for desired spreadFreqHz
+
+  spread_FsHz = resample(spread_lowFs(Ntaps+1:Nsam_low), FsHz, lowFs);
+
+  % return some states for optional unit testing
+  states.x = x;
+  states.y = y;
+  states.b = b;
+
+endfunction
+
+

diff --git a/codec2-dev/octave/doppler_spread_ut.m b/codec2-dev/octave/doppler_spread_ut.m
new file mode 100644 (file)
index 0000000..f7d96b0
--- /dev/null
+++ b/codec2-dev/octave/doppler_spread_ut.m
@@ -0,0 +1,51 @@
+% doppler_spread_ut.m
+% David Rowe Jan 2016
+%
+% Unit test script for doppler_spread
+
+f = 1;
+Fs = 8000;
+N  = Fs*10;
+
+[spread states] = doppler_spread(f, Fs, N);
+
+% use spreading samples to modulate 1000Hz sine wave
+% You can listen to this with: sine1k_1Hz.raw
+
+%   $ play -t raw -r 8000 -s -2 
+s = cos(2*pi*(1:N)*1000/Fs);
+s = s .* spread;
+s = real(s)*5000;
+fs = fopen("sine1k_1Hz.raw","wb"); fwrite(fs,s,"short"); fclose(fs);
+
+% Some plots
+
+x = states.x; y = states.y; b = states.b;
+
+H = freqz(b,1,x);
+
+figure(1)
+clf
+subplot(211)
+plot(x,y,';target;')
+title('Gaussian Filter Freq Resp Lin');
+legend('boxoff');
+subplot(212)
+plot(x,20*log10(y),';target;')
+hold on;
+plot(x,20*log10(y),'g+;actual;')
+hold off;
+axis([0 f*10/2 -60 0])
+title('Gaussian Filter Freq Resp dB');
+xlabel('Freq (Hz)');
+legend('boxoff');
+
+figure(2);
+subplot(211)
+plot(abs(spread))
+title('Spreading Function Magnitude');
+subplot(212)
+plot(s)
+title('1000Hz Sine Wave');
+xlabel('Time (samples)')
+