+++ /dev/null
-README.txt
-For codec2/octave directory
-Created 24 June 2012 by David Rowe
-
-1/ To support some of the Octave scripts (in particular fdmdv) in this
- directory the following Octave packages need to be installed:
-
- control image miscellaneous optim signal specfun struct
-
-2/ Download these tar balls from:
-
- http://octave.sourceforge.net/packages.php
-
-3/ Install each package from the Octave command line with:
-
- octave:3> pkg install package_file_name.tar.gz
-
end
XdB /= num_frames;
-XdB -= max(XdB);
+XdB -= max(20*log10(N));
figure(1)
+clf
plot((0:N/2-1)*Fs/(1000*N), XdB)
grid
ylabel('Amplitude dB')
xlabel('Frequency (kHz)');
-axis([0 Fs/(2*1000) -80 0])
+axis([0 Fs/(2*1000) -30 80])
+
+figure(2)
+clf
+plot(s)
+
% adcres.m
% David Rowe 18 Feb 2015
%
-% ADC resamping simulation
+% ADC resamping simulation, IIR tuner development.
% [ ] quantisation of ADC
% [ ] SNR at ADC input, SNR at resampler output
t = (0:(fs-1));
M = 45;
beta1 = 0.999;
-beta2 = 1 - (1-0.999)*M;
+beta2 = 1 - (1-beta1)*M;
s1 = [fs zeros(1,fs-1)]; % noise floor, continuous interferers
s2 = 100*4*cos(t*2*pi*f2/fs); % wanted signal 40dB above interferers
--- /dev/null
+% load_comp.m
+% David Rowe Sep 2015
+
+function s = load_comp(fn)
+ fs=fopen(fn,"rb");
+ s = fread(fs,Inf,"float32");
+ ls = length(s);
+ s = s(1:2:ls) + j*s(2:2:ls);
+endfunction
--- /dev/null
+% mag_to_phase.m
+%
+% David Rowe Sep 2015
+%
+% Slighly modified version of http://www.dsprelated.com/showcode/20.php
+%
+% Given a magnitude spectrum in dB, returns a minimum-phase phase
+% spectra. Both must be sampled at a Nfft. My understanding of this
+% is rather dim, but a working example is good place to start!
+
+
+function [phase s] = mag_to_phase(Gdbfk)
+ Nfft = 512; % FFT size to use
+
+ Ns = length(Gdbfk); if Ns~=Nfft/2+1, error("confusion"); end
+ Sdb = [Gdbfk,Gdbfk(Ns-1:-1:2)]; % install negative-frequencies
+
+ S = 10 .^ (Sdb/20); % convert to linear magnitude
+ s = ifft(S); % desired impulse response
+ s = real(s); % any imaginary part is quantization noise
+ tlerr = 100*norm(s(round(0.9*Ns:1.1*Ns)))/norm(s);
+ disp(sprintf([' Time-limitedness check: Outer 20%% of impulse ' ...
+ 'response is %0.2f %% of total rms'],tlerr));
+ % = 0.02 percent
+
+ if tlerr>1.0 % arbitrarily set 1% as the upper limit allowed
+ disp(' Increase Nfft and/or smooth Sdb\n');
+ end
+
+ c = ifft(Sdb); % compute real cepstrum from log magnitude spectrum
+
+ % Check aliasing of cepstrum (in theory there is always some):
+
+ caliaserr = 100*norm(c(round(Ns*0.9:Ns*1.1)))/norm(c);
+ disp(sprintf([' Cepstral time-aliasing check: Outer 20%% of ' ...
+ 'cepstrum holds %0.2f %% of total rms\n'],caliaserr));
+
+ if caliaserr>1.0 % arbitrary limit
+ disp(' Increase Nfft and/or smooth Sdb to shorten cepstrum\n');
+ end
+
+ % Fold cepstrum to reflect non-min-phase zeros inside unit circle:
+
+ cf = [c(1), c(2:Ns-1)+c(Nfft:-1:Ns+1), c(Ns), zeros(1,Nfft-Ns)];
+ Cf = fft(cf); % = dB_magnitude + j * minimum_phase
+
+ % The maths says we are meant to be using log(x), not 20*log10(x),
+ % so we need to sclae the phase toa ccount for this:
+ % log(x) = 20*log10(x)/scale;
+
+ scale = (20/log(10));
+ phase = imag(Cf)/scale;
+endfunction
+
Gdbfk = 20*log10(abs(Aw(1:Nfft/2+1)));
end
- Ns = length(Gdbfk); if Ns~=Nfft/2+1, error("confusion"); end
- Sdb = [Gdbfk,Gdbfk(Ns-1:-1:2)]; % install negative-frequencies
-
- S = 10 .^ (Sdb/20); % convert to linear magnitude
- s = ifft(S); % desired impulse response
- s = real(s); % any imaginary part is quantization noise
- tlerr = 100*norm(s(round(0.9*Ns:1.1*Ns)))/norm(s);
- disp(sprintf([' Time-limitedness check: Outer 20%% of impulse ' ...
- 'response is %0.2f %% of total rms'],tlerr));
- % = 0.02 percent
-
- if tlerr>1.0 % arbitrarily set 1% as the upper limit allowed
- disp(' Increase Nfft and/or smooth Sdb\n');
- end
-
- c = ifft(Sdb); % compute real cepstrum from log magnitude spectrum
-
- % Check aliasing of cepstrum (in theory there is always some):
+ [phase s] = mag_to_phase(Gdbfk);
- caliaserr = 100*norm(c(round(Ns*0.9:Ns*1.1)))/norm(c);
- disp(sprintf([' Cepstral time-aliasing check: Outer 20%% of ' ...
- 'cepstrum holds %0.2f %% of total rms\n'],caliaserr));
-
- if caliaserr>1.0 % arbitrary limit
- disp(' Increase Nfft and/or smooth Sdb to shorten cepstrum\n');
- end
-
- % Fold cepstrum to reflect non-min-phase zeros inside unit circle:
-
- cf = [c(1), c(2:Ns-1)+c(Nfft:-1:Ns+1), c(Ns), zeros(1,Nfft-Ns)];
- Cf = fft(cf); % = dB_magnitude + j * minimum_phase
-
- phase = imag(Cf)/(20/log(10));
endfunction
+
sim_dqpsk_hf = ber_test(sim_in, 'dqpsk');\r
sim_in.ldpc_code_rate = 1/2;\r
sim_in.ldpc_code = 1;\r
+\r
+ sim_in.hf_sim = 0;\r
+ sim_qpsk_ldpc = ber_test(sim_in, 'qpsk');\r
+\r
+ sim_in.hf_sim = 1;\r
sim_qpsk_hf_ldpc = ber_test(sim_in, 'qpsk');\r
sim_in.hf_mag_only = 0;\r
sim_dqpsk_hf_ldpc = ber_test(sim_in, 'dqpsk');\r
hold on;\r
semilogy(sim_qpsk.Ebvec, sim_qpsk.BERvec,'g;QPSK AWGN;')\r
semilogy(sim_qpsk_hf.Ebvec, sim_qpsk_hf.BERvec,'r;QPSK HF;')\r
+ semilogy(sim_qpsk_ldpc.Ebvec, sim_qpsk_ldpc.BERldpcvec,'bk;QPSK HF;')\r
semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'g;DQPSK AWGN;')\r
semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'r;DQPSK HF;')\r
semilogy(sim_qpsk_hf_ldpc.Ebvec, sim_qpsk_hf_ldpc.BERldpcvec,'b;QPSK HF LDPC 1/2;')\r
demod_type = 0;
decoder_type = 0;
max_iterations = 100;
- EsNo = 10; % decoder needs an estimated channel EsNo (linear ratio, not dB)
code_param = ldpc_init(aqpsk.code_rate, aqpsk.framesize, modulation, mod_order, mapping);
assert(code_param.data_bits_per_frame == aqpsk.Ndatabits);
clf
semilogy(EbNodB, uncoded_BER_theory,'r;uncoded QPSK theory;')
hold on;
-semilogy(EbNodB, sim_out.BER_uc,'r+;uncoded QPSK actual;')
-semilogy(EbNodB-10*log10(sim_out.code_rate), sim_out.BER+1E-10,'g;LDPC coded QPSK simulation;');
+semilogy(EbNodB, sim_out.BER_uc,'r+;uncoded QPSK simulated;')
+semilogy(EbNodB-10*log10(sim_out.code_rate), sim_out.BER+1E-10,'g;LDPC coded QPSK simulated;');
hold off;
grid('minor')
xlabel('Eb/No (dB)')
ylabel('BER')
axis([min(EbNodB) max(EbNodB) min(uncoded_BER_theory) 1])
-% add noise
% correct freq offset
% output "demodulated" tx frame to a file of bits for 3rd party modulator
-% actually estimate Es/No
-% meas impl loss for pilot symbols
-% put in Es/No loop, measure coded and uncoded errors
-% separate into end and dec
-
+% actually estimate Es/No for LDPC
+% improve phase est impl loss
+% separate into enc and demod functions
+% generate filter coeffs, save to C header file
projects / freetel-svn-tracking.git / commitdiff