From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Thu, 17 Sep 2015 20:47:35 +0000 (+0000)
Subject: stand alone magnitude to phase function
X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=ac5a0d8fee0fbdb291d75764229f7185be73ba1c;p=freetel-svn-tracking.git

stand alone magnitude to phase function

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

diff --git a/codec2-dev/octave/README.txt b/codec2-dev/octave/README.txt
deleted file mode 100644
index 05d63f20..00000000
--- a/codec2-dev/octave/README.txt
+++ /dev/null
@@ -1,17 +0,0 @@
-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
-
diff --git a/codec2-dev/octave/adc_sfdr_ut.m b/codec2-dev/octave/adc_sfdr_ut.m
index c2363ed9..06fcf307 100644
--- a/codec2-dev/octave/adc_sfdr_ut.m
+++ b/codec2-dev/octave/adc_sfdr_ut.m
@@ -17,12 +17,18 @@ for i=1:num_frames
 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)
+
 
diff --git a/codec2-dev/octave/adcres.m b/codec2-dev/octave/adcres.m
index f4f04ac0..69dae0d4 100644
--- a/codec2-dev/octave/adcres.m
+++ b/codec2-dev/octave/adcres.m
@@ -1,7 +1,7 @@
 % 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
@@ -19,7 +19,7 @@ f4 = f1 - 207E3;
 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
diff --git a/codec2-dev/octave/load_comp.m b/codec2-dev/octave/load_comp.m
new file mode 100644
index 00000000..b9fa6866
--- /dev/null
+++ b/codec2-dev/octave/load_comp.m
@@ -0,0 +1,9 @@
+% 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
diff --git a/codec2-dev/octave/mag_to_phase.m b/codec2-dev/octave/mag_to_phase.m
new file mode 100644
index 00000000..4e553404
--- /dev/null
+++ b/codec2-dev/octave/mag_to_phase.m
@@ -0,0 +1,54 @@
+% 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
+
diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m
index 80dc173d..b170810f 100644
--- a/codec2-dev/octave/newamp.m
+++ b/codec2-dev/octave/newamp.m
@@ -477,38 +477,8 @@ function [phase Gdbfk s Aw] = determine_phase(model, f, ak)
     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
 
+
diff --git a/codec2-dev/octave/test_qpsk2.m b/codec2-dev/octave/test_qpsk2.m
index bf7bd8e9..dbb43cbf 100644
--- a/codec2-dev/octave/test_qpsk2.m
+++ b/codec2-dev/octave/test_qpsk2.m
@@ -527,6 +527,11 @@ function ideal
   sim_dqpsk_hf            = ber_test(sim_in, 'dqpsk');
   sim_in.ldpc_code_rate   = 1/2;
   sim_in.ldpc_code        = 1;
+
+  sim_in.hf_sim           = 0;
+  sim_qpsk_ldpc           = ber_test(sim_in, 'qpsk');
+
+  sim_in.hf_sim           = 1;
   sim_qpsk_hf_ldpc        = ber_test(sim_in, 'qpsk');
   sim_in.hf_mag_only      = 0;
   sim_dqpsk_hf_ldpc       = ber_test(sim_in, 'dqpsk');
@@ -537,6 +542,7 @@ function ideal
   hold on;
   semilogy(sim_qpsk.Ebvec, sim_qpsk.BERvec,'g;QPSK AWGN;')
   semilogy(sim_qpsk_hf.Ebvec, sim_qpsk_hf.BERvec,'r;QPSK HF;')
+  semilogy(sim_qpsk_ldpc.Ebvec, sim_qpsk_ldpc.BERldpcvec,'bk;QPSK HF;')
   semilogy(sim_dqpsk.Ebvec, sim_dqpsk.BERvec,'g;DQPSK AWGN;')
   semilogy(sim_dqpsk_hf.Ebvec, sim_dqpsk_hf.BERvec,'r;DQPSK HF;')
   semilogy(sim_qpsk_hf_ldpc.Ebvec, sim_qpsk_hf_ldpc.BERldpcvec,'b;QPSK HF LDPC 1/2;')
diff --git a/codec2-dev/octave/tqpsk.m b/codec2-dev/octave/tqpsk.m
index 099a454e..4ec8765a 100644
--- a/codec2-dev/octave/tqpsk.m
+++ b/codec2-dev/octave/tqpsk.m
@@ -62,7 +62,6 @@ function sim_out = run_simulation(sim_in)
   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);
@@ -151,19 +150,17 @@ figure(1)
 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