From: drowe67 Date: Tue, 6 Jun 2017 04:04:11 +0000 (+0000) Subject: Attempt at correcting alias errors in rate K resampling, sig improvement for hst1a X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=4809a8a99b1b507ef50e44c26f79c6c7092fedc9;p=freetel-svn-tracking.git Attempt at correcting alias errors in rate K resampling, sig improvement for hst1a git-svn-id: https://svn.code.sf.net/p/freetel/code@3157 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/newamp.m b/codec2-dev/octave/newamp.m index 01089342..a1c64bc2 100644 --- a/codec2-dev/octave/newamp.m +++ b/codec2-dev/octave/newamp.m @@ -46,6 +46,43 @@ function y = interp_para(xp, yp, x) endfunction +% choose largest sample in band, idea is we care more about finding +% peaks, can handle some error in frequency. x are non linear +% (arbitrary) sampling points in kHz + +function y = interp_largest(f0_Hz, AmdB, x_kHz) + L = length(AmdB); + x = x_kHz*1000; + y = zeros(1,length(x)); + bw = x(2) - x(1); + k = 1; + + for i=1:length(x) + + % determine limits of this band + + if i>1 + bw = x(i) - x(i-1); + end + band_low = x(i) - bw/2; band_high = x(i) + bw/2; + + % map band limits to harmonics + + if x(i) < f0_Hz + m_low = m_high = 1; + else + m_low = round(band_low/f0_Hz); m_high = round(band_high/f0_Hz)-1; + m_low = max(1, m_low); m_high = min(L, m_high); m_high = max(m_low, m_high); + end + + printf("L: %d f0: %f i: %d band_low: %f band_high: %f m_low: %d m_high: %d\n",L, f0_Hz, i, band_low, band_high, m_low, m_high); + % find max in band + + y(i) = max(AmdB(m_low:m_high)); + end + +endfunction + % simple linear interpolator function y = interp_linear(xp, yp, x) @@ -172,7 +209,7 @@ endfunction function rate_K_sample_freqs_kHz = mel_sample_freqs_kHz(K) - mel_start = ftomel(200); mel_end = ftomel(3700); + mel_start = ftomel(100); mel_end = ftomel(0.95*4000); step = (mel_end-mel_start)/(K-1); mel = mel_start:step:mel_end; rate_K_sample_freqs_Hz = 700*((10 .^ (mel/2595)) - 1); @@ -222,6 +259,64 @@ function [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f(model, endfunction +function [rate_K_vec_corrected orig_error error nasty_error_log nasty_error_m_log] = correct_rate_K_vec(rate_K_vec, rate_K_sample_freqs_kHz, AmdB, AmdB_, K, Wo, L, Fs) + + % aliasing correction -------------------------------------- + + % The mel sample rate decreases as frequency increases. Look for + % any regions above 1000Hz where we have missed definition of a + % spectral peak (formant) due to aliasing. Adjust the rate K + % sample levels to restore peaks. Theory is that correct + % definition of a formant is less important than the frequency of + % the formant. As long as we define a formant in that general + % frequency area it will sound OK. + + Am_freqs_kHz = (1:L)*Wo*4/pi; + + % Lets see where we have made an error + + error = orig_error = AmdB - AmdB_; + + Ncorrections = 3; % maximum number of rate K samples to correct + error_thresh = 3; % only worry about errors larger than thresh + + start_m = floor(L*1000/(Fs/2)); + error(1:start_m) = 0; % first 1000Hz is densly sampled so ignore + nasty_error_m_log = []; nasty_error_log = []; + + + rate_K_vec_corrected = rate_K_vec; + for i=1:Ncorrections + [mx mx_m] = max(error); + + if mx > error_thresh + nasty_error_log = [nasty_error_log mx]; + nasty_error_m_log = [nasty_error_m_log mx_m]; + + % find closest rate K sample to nasty error + + nasty_error_freq = mx_m*Wo*Fs/(2*pi*1000) + [tmp closest_k] = min(abs(rate_K_sample_freqs_kHz - nasty_error_freq)); + rate_K_vec_corrected(closest_k) = AmdB(mx_m); + + % zero out error in this region and look for another large error region + + k = max(1, closest_k-1); + rate_K_prev_sample_kHz = rate_K_sample_freqs_kHz(k) + k = min(K, closest_k+1); + rate_K_next_sample_kHz = rate_K_sample_freqs_kHz(k) + + [tmp st_m] = min(abs(Am_freqs_kHz - rate_K_prev_sample_kHz)); + [tmp en_m] = min(abs(Am_freqs_kHz - rate_K_next_sample_kHz)); + if closest_k == K + en_m = L; + end + error(st_m:en_m) = 0; + end + end +endfunction + + % Take a rate K surface and convert back to time varying rate L function [model_ AmdB_] = resample_rate_L(model, rate_K_surface, rate_K_sample_freqs_kHz) @@ -1266,3 +1361,12 @@ function lmin = abys(AmdB_, AmdB, Wo, L, mask_sample_freqs_kHz) endfunction +function rate_K_surface_no_slope = remove_slope(rate_K_surface) + [frames K] = size(rate_K_surface); + rate_K_surface_no_slope = zeros(frames,K); + for f=1:frames + [gradient intercept] = linreg(1:K, rate_K_surface(f,:), K); + printf("f: %d gradient: %f intercept: %f\n", f, gradient, intercept); + rate_K_surface_no_slope(f,:) = rate_K_surface(f,:) - (intercept + gradient*(1:K)); + end +endfunction diff --git a/codec2-dev/octave/newamp1_batch.m b/codec2-dev/octave/newamp1_batch.m index 58793139..3dbbbc7c 100644 --- a/codec2-dev/octave/newamp1_batch.m +++ b/codec2-dev/octave/newamp1_batch.m @@ -7,7 +7,7 @@ #{ Octave script to batch process model parameters using the new - amplitude model, version 1. Outputs anotehr set of model parameters + amplitude model, version 1. Outputs another set of model parameters that can be fed to c2sim for listening tests. The companion newamp1_fbf.m script is used to visualise the processing frame by frame @@ -44,7 +44,7 @@ function surface = newamp1_batch(input_prefix, output_prefix) newamp; more off; - max_amp = 80; + max_amp = 160; postfilter = 0; % optional postfiler that runs on Am, not used atm synth_phase = 1; @@ -68,12 +68,16 @@ function surface = newamp1_batch(input_prefix, output_prefix) %model_ = experiment_filter(model); %model_ = experiment_filter_dec_filter(model); - %[model_ surface] = experiment_mel_freq(model, 1, 1, voicing); + [model_ surface] = experiment_mel_freq(model, 0, 1, voicing); + %model_ = experiment_dec_abys(model, 8, 1, 1, 1, voicing); +#{ [model_ voicing_ indexes] = experiment_rate_K_dec(model, voicing); % encoder/decoder, lets toss away results except for indexes - %[model_ voicing_] = model_from_indexes(indexes); % decoder uses just indexes, outputs vecs for synthesis [model_ voicing_] = model_from_indexes_fbf(indexes); % decoder uses just indexes, outputs vecs for synthesis +#} + + %[model_ voicing_] = model_from_indexes(indexes); % decoder uses just indexes, outputs vecs for synthesis %model_ = experiment_dec_linear(model_); %model_ = experiment_energy_rate_linear(model, 1, 0); @@ -121,7 +125,7 @@ function surface = newamp1_batch(input_prefix, output_prefix) % synthesis phase spectra from magnitiude spectra using minimum phase techniques - fft_enc = 128; + fft_enc = 256; phase = determine_phase(model_, f, fft_enc); assert(length(phase) == fft_enc); %Aw = zeros(1, fft_enc*2); @@ -492,9 +496,21 @@ endfunction function [model_ rate_K_surface] = experiment_mel_freq(model, vq_en=0, plots=1, voicing) [frames nc] = size(model); - K = 20; - [rate_K_surface rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model, K); - + K = 20; Fs = 8000; + + for f=1:frames + Wo = model(f,1); + L = model(f,2); + Am = model(f,3:(L+2)); + AmdB = 20*log10(Am); + Am_freqs_kHz = (1:L)*Wo*4/pi; + [rate_K_vec rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model(f,:), K); + [tmp_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec, rate_K_sample_freqs_kHz); + [rate_K_vec_corrected orig_error error nasty_error_log nasty_error_m_log] = correct_rate_K_vec(rate_K_vec, rate_K_sample_freqs_kHz, AmdB, AmdB_, K, Wo, L, Fs); + %rate_K_surface(f,:) = rate_K_vec; + rate_K_surface(f,:) = rate_K_vec_corrected; + end + if plots figure(1); clf; mesh(rate_K_surface); figure(2); clf; plot_dft_surface(rate_K_surface) @@ -1057,3 +1073,4 @@ function plot_dft_surface(surface) Fs = 100; dF = Fs/frames; mesh(1:K, (1:frames/2)*dF, 20*log10(abs(dft_surface(1:frames/2,:,:)))); endfunction + diff --git a/codec2-dev/octave/newamp1_fbf.m b/codec2-dev/octave/newamp1_fbf.m index 2946d5de..fe7bb2ba 100644 --- a/codec2-dev/octave/newamp1_fbf.m +++ b/codec2-dev/octave/newamp1_fbf.m @@ -18,10 +18,10 @@ function newamp1_fbf(samname, f=10) newamp; more off; - quant_en = 0; pf_en = 0; plot_phase = 1; + quant_en = 0; pf_en = 0; plot_phase = 0; melvq; - - K=20; load train_120_vq; m=5; + + Fs = 8000; K=20; load train_120_vq; m=5; % load up text files dumped from c2sim --------------------------------------- @@ -66,9 +66,14 @@ function newamp1_fbf(samname, f=10) AmdB = 20*log10(Am); Am_freqs_kHz = (1:L)*Wo*4/pi; - % plots for mel sampling + % resample at rate K [rate_K_vec rate_K_sample_freqs_kHz] = resample_const_rate_f_mel(model(f,:), K); + [tmp_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec, rate_K_sample_freqs_kHz); + [rate_K_vec_corrected orig_error error nasty_error_log nasty_error_m_log] = correct_rate_K_vec(rate_K_vec, rate_K_sample_freqs_kHz, AmdB, AmdB_, K, Wo, L, Fs); + [tmp_ AmdB_corrected] = resample_rate_L(model(f,:), rate_K_vec, rate_K_sample_freqs_kHz); + + % quantisation simulation --------------------------------- mean_f = mean(rate_K_vec); rate_K_vec_no_mean = rate_K_vec - mean_f; @@ -85,6 +90,7 @@ function newamp1_fbf(samname, f=10) rate_K_vec_ = rate_K_vec_no_mean_ + mean_f; [model_ AmdB_] = resample_rate_L(model(f,:), rate_K_vec_, rate_K_sample_freqs_kHz); + %[model_ AmdB_largest] = resample_rate_L(model(f,:), rate_K_largest, rate_K_sample_freqs_kHz); % plots ---------------------------------- @@ -94,17 +100,20 @@ function newamp1_fbf(samname, f=10) axis([1 4000 -20 80]); hold on; + plot((1:L)*Wo*4000/pi, AmdB,";Am;b+-"); - plot(rate_K_sample_freqs_kHz*1000, rate_K_vec, ';rate K mel;g+-'); - if quant_en >= 1 - plot((1:L)*Wo*4000/pi, AmdB_,";Am quant;k+-"); - end - if quant_en == 2 - plot(rate_K_sample_freqs_kHz*1000, rate_K_vec_, ';rate K mel quant;r+-'); - end + stem(rate_K_sample_freqs_kHz*1000, rate_K_vec, 'g'); + plot(rate_K_sample_freqs_kHz*1000, rate_K_vec, ';rate K;g:'); + plot(rate_K_sample_freqs_kHz*1000, rate_K_vec_corrected, ';rate K;g-', 'linewidth', 2); + + plot((1:L)*Wo*4000/pi, -10 + orig_error ,";AmdB - AmdB_ error;k+-"); + plot((1:L)*Wo*4000/pi, -10 + error ,";AmdB - AmdB_ error;r-"); + plot(nasty_error_m_log*Wo*4000/pi, AmdB(nasty_error_m_log), 'ro', 'markersize', 10); + hold off; + #{ figure(3); clf; title('Frequency Domain 2'); @@ -116,6 +125,7 @@ function newamp1_fbf(samname, f=10) plot(rate_K_sample_freqs_kHz*1000, rate_K_vec_no_mean_, ';rate K mel no mean quant;r+-'); end hold off; + #} if plot_phase phase_512 = determine_phase(model_, 1, 512); @@ -135,7 +145,7 @@ function newamp1_fbf(samname, f=10) if k == 'm' quant_en++; - if quant_en > 2 + if quant_en > 1 quant_en = 0; end endif