From: drowe67 Date: Thu, 17 Sep 2015 12:44:40 +0000 (+0000) Subject: generating phase spectra from mask, works quite well on hts1a, hts2a, ve9qrp_10s X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=1b7a817de3f761b5681b88a5e9071ba1a67841cd;p=freetel-svn-tracking.git 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 bbe72f50..80dc173d 100644 --- 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 50ca81cb..591addc9 100644 --- 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 4a90df8a..119a5814 100644 --- 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