From: drowe67 Date: Mon, 25 Jan 2016 00:29:49 +0000 (+0000) Subject: doppler spreading functions for HF channel simulation X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=36e118f4e9108a5d9e89ad184f9184e1f3d8fd60;p=freetel-svn-tracking.git 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 index 00000000..4bc6c62a --- /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 index 00000000..f7d96b01 --- /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)') +