generating phase spectra from mask, works quite well on hts1a, hts2a, ve9qrp_10s

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

diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m
index bbe72f5062a7185a223a732e6e1ae5a2b5acf151..80dc173dd219f3bf9a150c9e21fb4c7c8a02822f 100644 (file)
--- a/codec2-dev/octave/newamp.m
+++ b/codec2-dev/octave/newamp.m
@@ -330,6 +330,9 @@ function maskdB = resonator(freq_tone_kHz, mask_sample_freqs_kHz)
   end
 endfunction
 
+
+% sampling the mask in one frame using an arbitrary set of samplng frequencies
+
 function maskdB = resample_mask(model, f, mask_sample_freqs_kHz)
   max_amp = 80;
 
@@ -342,7 +345,7 @@ function maskdB = resample_mask(model, f, mask_sample_freqs_kHz)
 endfunction
 
 
-% decimate frame rate of mask, use interpolation in the log domain 
+% decimate frame rate of mask, use linear interpolation in the log domain 
 
 function maskdB_ = decimate_frame_rate(maskdB, model, decimate, f, frames, mask_sample_freqs_kHz)
     max_amp = 80;
@@ -449,4 +452,63 @@ function amp_scatter(samname)
   hist(ampsdB,20)
 endfunction
 
-%amp_scatter("../build_linux/src/all")
+
+% Determine a phase spectra from a magnitude spectra
+% from http://www.dsprelated.com/showcode/20.php
+% Haven't _quite_ figured out how this works but have to start somewhere ....
+%
+% TODO: we may be able to sample at a lower rate, like mWo
+%       but start with something that works
+
+function [phase Gdbfk s Aw] = determine_phase(model, f, ak)
+  Nfft    = 512;  % FFT size to use 
+  Fs      = 8000;
+  max_amp = 80;
+  L       = min([model(f,2) max_amp-1]);
+  Wo      = model(f,1);
+
+  mask_sample_freqs_kHz = (Fs/1000)*[0:Nfft/2]/Nfft;           % fft frequency grid (nonneg freqs)
+  Gdbfk = resample_mask(model, f, mask_sample_freqs_kHz);
+
+  % optional input of aks for testing
+
+  if nargin == 3
+    Aw = 1 ./ fft(ak,Nfft);
+    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):
+
+  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/newamp_batch.m b/codec2-dev/octave/newamp_batch.m
index 50ca81cbd817fa9aa67018bac76a7a913717adc5..591addc9c03715bd7451317196ec4cb6b7f1b639 100644 (file)
--- a/codec2-dev/octave/newamp_batch.m
+++ b/codec2-dev/octave/newamp_batch.m
@@ -16,7 +16,7 @@
 
 % process a whole file and write results
 
-function newamp_batch(samname, optional_Am_out_name)
+function newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name)
   newamp;
   more off;
 
@@ -31,12 +31,17 @@ function newamp_batch(samname, optional_Am_out_name)
 
   if nargin == 1
     Am_out_name = sprintf("%s_am.out", samname);
+    Aw_out_name = sprintf("%s_aw.out", samname);
   end
-  if nargin == 2
+  if nargin >= 2
     Am_out_name = optional_Am_out_name;
   end
-    
-  fam = fopen(Am_out_name,"wb"); 
+  if nargin >= 3
+    Aw_out_name = optional_Aw_out_name;
+  end
+   
+  fam  = fopen(Am_out_name,"wb"); 
+  faw = fopen(Aw_out_name,"wb"); 
 
   for f=1:frames
     printf("frame: %d", f);
@@ -87,9 +92,19 @@ function newamp_batch(samname, optional_Am_out_name)
     Am_ = zeros(1,max_amp);
     Am_(2:L) = 10 .^ (maskdB_pf(1:L-1)/20);  % C array doesnt use A[0]
     fwrite(fam, Am_, "float32");
+
+    fft_enc = 512;
+    phase = determine_phase(model, f);
+    assert(length(phase) == fft_enc);
+    Aw = zeros(1, fft_enc*2); 
+    Aw(1:2:fft_enc*2) = cos(phase);
+    Aw(2:2:fft_enc*2) = -sin(phase);
+ 
+    fwrite(faw, Aw, "float32");    
   end
 
   fclose(fam);
+  fclose(faw);
 
   printf("\n");
 

diff --git a/codec2-dev/octave/newamp_fbf.m b/codec2-dev/octave/newamp_fbf.m
index 4a90df8a0c27409ce9463b61c2a2e4e732968e4e..119a5814639759942ae84a144fc9651a3257165d 100644 (file)
--- a/codec2-dev/octave/newamp_fbf.m
+++ b/codec2-dev/octave/newamp_fbf.m
@@ -93,8 +93,8 @@ function newamp_fbf(samname, f)
 
     % decimated in time
 
-    maskdB = decimate_frame_rate(maskdB, model, 4, f, frames, mask_sample_freqs_kHz);
-    plot(mask_sample_freqs_kHz*1000, maskdB, 'k');
+    %maskdB = decimate_frame_rate(maskdB, model, 4, f, frames, mask_sample_freqs_kHz);
+    %plot(mask_sample_freqs_kHz*1000, maskdB, 'k');
 
     % optionally plot all masking curves
 
@@ -108,14 +108,23 @@ function newamp_fbf(samname, f)
 
     hold off;
 
+    ak_name = strcat(samname,"_ak.txt");
+    ak = load(ak_name);
+
+    [phase Sdb s Aw] = determine_phase(model, f, ak(f,:));
     figure(3)
-    plot(target_log,'g')
+    subplot(211)
+    plot(Sdb)
+    hold on;
+    plot(20*log10(abs(Aw(1:256))),'g')
+    hold off;
+    subplot(212)
+    plot(phase(1:256))
     hold on;
-    plot(candidate_log(3,:),'b');
-    plot(candidate_log(5,:),'b');
-    plot(error_log(3,:),'r');
-    plot(error_log(5,:),'r');
+    plot(angle(Aw(1:256)),'g')
     hold off;
+    figure(4)
+    plot(s)
 
     % interactive menu