1;
melvq;
+
+function [maskdB_ maskdB_cyclic Dabs dk_] = decimate_in_freq(maskdB, cyclic=1, Nkeep=7)
+
+ % Lets try to come up with a smoothed, cyclic model. Replace
+ % points from 3500 Hz to 4000Hz with a sequence that joins without
+ % a step to points at the 0Hz end of the spectrum. This will make
+ % it more cyclical and make the DFT happier, less high freq
+ % energy. Yes, happier is an extremely technical term.
+
+ L = length(maskdB);
+ anchor = floor(7*L/8);
+ xpts = [ anchor-1 anchor L+1 L+2];
+ ypts = [ maskdB(anchor-1) maskdB(anchor) maskdB(1) maskdB(2)];
+ mask_pp = splinefit(xpts, ypts, 1);
+ maskdB_cyclic = [maskdB(1:anchor) ppval(mask_pp, anchor+1:L)];
+
+ % Now DFT, truncating DFT coeffs to undersample
+
+ if cyclic
+ D = fft(maskdB_cyclic);
+ else
+ D = fft(maskdB);
+ end
+ Dabs = abs(D); % this returned for plotting
+ D_ = D; % + 10*randn(1,L) + 10*j*randn(1,L);
+ D_(Nkeep+1:L-Nkeep) = 0; % truncate
+ d = ifft(D_); % back to spectrum at rate L
+ maskdB_ = real(d);
+ Dk_ = [0 D(2:Nkeep-1) real(D(Nkeep)) D(L-Nkeep+1:L)]; % build rate Nkeep vector for quantisation
+ dk_ = real(ifft(Dk_));
+
+ % OK now fix up last 500Hz, taper down 10dB at 4000Hz
+
+ xpts = [ anchor-1 anchor L];
+ ypts = [ maskdB_(anchor-1) maskdB_(anchor) (maskdB_(anchor)-10)];
+ mask_pp = splinefit(xpts, ypts, 1);
+ maskdB_ = [maskdB_(1:anchor) ppval(mask_pp, anchor+1:L)];
+
+endfunction
+
+
% Create a "decimated mask" model using just a few samples of
% critical band filter centre frequencies. For voiced speech,
% we fit the amplitude of these samples to a straight line.
% process a whole file and write results
-function [non_masked_f_log non_masked_amp_log] = newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name)
+function dk_log = newamp_batch(samname, optional_Am_out_name, optional_Aw_out_name)
newamp;
more off;
max_amp = 80;
- decimate_in_freq = 1;
+ dec_in_freq = 1;
postfilter = 1;
- decimate_in_time = 0;
- synth_phase = 1;
+ dec_in_time = 0;
+ synth_phase = 0;
freq_quant = 0;
amp_quant = 0;
- non_masked_f_log = [];
- non_masked_m_log = [];
- non_masked_amp_log = [];
- best_min_mse_log = [];
+ dk_log = [];
model_name = strcat(samname,"_model.txt");
model = load(model_name);
%pp_bw = gen_pp_bw;
+ sd_sum = 0;
for f=1:frames
printf("%d ", f);
model_(f,2) = L = min([model(f,2) max_amp-1]);
mask_sample_freqs_kHz = (1:L)*Wo*4/pi;
maskdB = mask_model(AmdB, Wo, L);
- if decimate_in_freq
- % decimate mask samples in frequency
- [decmaskdB masker_freqs_kHz masker_amps_dB min_error mse_log1 best_min_mse ind_log] = make_decmask_mel(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz, freq_quant, amp_quant);
- best_min_mse_log = [best_min_mse_log best_min_mse];
- ind_log
+ % save post filter positions for decode
- non_masked_amp = masker_amps_dB ;
- non_masked_amp_log = [non_masked_amp_log; non_masked_amp];
+ maskdB_ = maskdB;
+ a_non_masked_m = find(AmdB > maskdB);
+ non_masked_m(f,1:length(a_non_masked_m)) = a_non_masked_m;
- % Save this for decoder, so it knows where to apply post filter
- % Basically the frequencies of the AbyS samples
-
- non_masked_m(f,:) = min(round(masker_freqs_kHz/(Wo*4/pi)),L);
- non_masked_m_log = [non_masked_m_log; non_masked_m(f,:)];
-
- non_masked_f_log = [non_masked_f_log; masker_freqs_kHz];
- maskdB_ = decmaskdB;
-
- else
- % just approximate decimation in freq by using those mask samples we can hear
- maskdB_ = maskdB;
- a_non_masked_m = find(AmdB > maskdB);
- non_masked_m(f,1:length(a_non_masked_m)) = a_non_masked_m;
+ if dec_in_freq
+ [maskdB_ tmp1 D dk_] = decimate_in_freq(maskdB, 1);
+ dk_log = [dk_log; dk_];
end
+ sd_sum += sum(maskdB - maskdB_);
Am_ = zeros(1,max_amp);
Am_ = 10 .^ (maskdB_(1:L-1)/20);
model_(f,3:(L+1)) = Am_;
end
+ printf("\nsd_sum: %5.2f\n", sd_sum/frames);
% decoder loop -----------------------------------------------------
maskdB_ = 20*log10(Am_);
mask_sample_freqs_kHz = (1:L)*Wo*4/pi;
- %maskdB_ = mask_model(maskdB_, Wo, L);
- if decimate_in_time
+ if dec_in_time
% decimate mask samples in time
maskdB_ = decimate_frame_rate(model_, 4, f, frames, mask_sample_freqs_kHz);
m = max(find(non_masked_m(f,:) > 0));
a_non_masked_m = non_masked_m(f,1:m);
-
end
% post filter - bump up samples by 6dB, reduce mask by same level to normalise gain
Aw(1:2:fft_enc*2) = cos(phase);
Aw(2:2:fft_enc*2) = -sin(phase);
- if 0
- % optional plotting to ensure amd and phase aligned on a few frames
- figure(1)
- clf;
- subplot(211)
- plot(mask_sample_freqs_kHz, maskdB_);
- hold on;
- plot(mask_sample_freqs_kHz, maskdB_pf,'g');
- hold off;
- axis([0 4 0 70])
-
- subplot(212)
- plot((0:(fft_enc/2)-1)*8000/fft_enc, phase(1:fft_enc/2))
- axis([0 4000 -pi pi])
- k = kbhit();
- end
fwrite(faw, Aw, "float32");
end
end
printf("\n");
- printf("mse std: %5.2f\n", std(best_min_mse_log));
-
endfunction
% $ cd ~/codec2-dev/octave
% octave:14> newamp_fbf("../build_linux/src/hts1a",50)
+
+
function newamp_fbf(samname, f=10)
-
- more off;
newamp;
- phase_stuff = 0;
- plot_not_masked = 0;
+ more off;
plot_spectrum = 1;
- freq_quant = 0;
- amp_quant = 0;
- decimate_in_freq = 0;
+
+ % load up text files dumped from c2sim ---------------------------------------
sn_name = strcat(samname,"_sn.txt");
Sn = load(sn_name);
-
sw_name = strcat(samname,"_sw.txt");
Sw = load(sw_name);
-
model_name = strcat(samname,"_model.txt");
model = load(model_name);
[frames tmp] = size(model);
- if phase_stuff
- ak_name = strcat(samname,"_ak.txt");
- ak = load(ak_name);
- end
-
- % pp_bw = gen_pp_bw;
+ % Keyboard loop --------------------------------------------------------------
- plot_all_masks = 0;
k = ' ';
do
figure(1);
Am = model(f,3:(L+2));
AmdB = 20*log10(Am);
- % plotting
-
axis([1 4000 -20 80]);
hold on;
if plot_spectrum
plot((1:L)*Wo*4000/pi, AmdB,";Am;r");
plot((1:L)*Wo*4000/pi, AmdB,"r+");
- %plot((0:255)*4000/256, Sw(f,:),";Sw;");
end
[maskdB Am_freqs_kHz] = mask_model(AmdB, Wo, L);
+ [tmp1 tmp2 D] = decimate_in_freq(maskdB, 0);
+ [maskdB_ maskdB_cyclic D_cyclic dk_] = decimate_in_freq(maskdB, 1);
- plot(Am_freqs_kHz*1000, maskdB, 'g');
-
- % Lets try to come up with a smoothed model. Replace points from 3500 Hz on with
- % those than match low end of spectrum. This will make it more cyclical and make
- % DFT happier, less high freq energy. Need a better way to describe that.
-
- anchor = floor(7*L/8);
- xpts = [ anchor-1 anchor L+1 L+2];
- ypts = [ maskdB(anchor-1) maskdB(anchor) maskdB(1) maskdB(2)];
- mask_pp = splinefit(xpts, ypts, 1);
- maskdB_smooth = [maskdB(1:anchor) ppval(mask_pp, anchor+1:L)];
-
- plot(Am_freqs_kHz*1000, maskdB_smooth, 'b');
-
- if decimate_in_freq
- % decimate in frequency
-
- mask_sample_freqs_kHz = (1:L)*Wo*4/pi;
- [decmaskdB masker_freqs_kHz masker_amps_dB min_error mse_log1] = make_decmask_mel(maskdB, AmdB, Wo, L, mask_sample_freqs_kHz, freq_quant, amp_quant);
-
- % find turning points - prototype for finding PF freqs when we decimate in time
-
- if 0
- d = decmaskdB(2:L) - decmaskdB(1:L-1);
- size(d)
- tp = [];
- for m=1:L-1
- if (d(m) > 0) && (d(m+1) < 0)
- tp = [tp m+1];
- end
- end
- end
-
- figure(2)
-
- tonef_kHz = masker_freqs_kHz;
- nlm = length(tonef_kHz);
-
- plot(tonef_kHz*1000, zeros(1,nlm), ';AbyS Mask Freqs;bk+');
- if decimate_in_freq
- plot(mask_sample_freqs_kHz*1000, decmaskdB, ';AbyS Mask;m');
- end
- %plot(tp*Wo*4000/pi, 5*ones(1,length(tp)), ';PF Freqs;r+');
- legend('boxoff');
- %plot(mask_sample_freqs_kHz*1000, min_error);
-
- end
-
- hold off;
+ plot(Am_freqs_kHz*1000, maskdB, ';mask;g');
+ plot(Am_freqs_kHz*1000, maskdB_cyclic, ';mask cyclic;b');
+ plot(Am_freqs_kHz*1000, maskdB_, ';mask trunated;c');
% lets get a feel for the "spectrum" of the smoothed spectral envelope
% this will give us a feel for how hard it is to code, ideally we would like
figure(3)
clf
- D = abs(fft(maskdB));
- en = floor(length(D)/2+1);
+ en = L/2+1;
stem(D(2:en),'g')
hold on;
- D_smooth = abs(fft(maskdB_smooth));
- stem(D_smooth(2:en),'b')
+ stem(D_cyclic(2:en),'b')
hold off;
% let plot the cumulative amount of energy in each DFT
figure(4)
clf
- plot(cumsum(D(2:en)/sum(D(2:en))),'g');
+ plot(cumsum(D(2:en)/sum(D(2:en))),';cumsum;g');
hold on;
- plot(cumsum(D_smooth(2:en)/sum(D_smooth(2:en))),'b');
+ plot(cumsum(D_cyclic(2:en)/sum(D_cyclic(2:en))),';cumsum cyclic;b');
hold off;
axis([1 L 0 1])
- % try truncating DFT coeffs and see what that looks like
-
- D = fft(maskdB_smooth);
- keep = 7;
- D(keep+1:L-keep) = 0;
- d = ifft(D);
-
- figure(2)
- hold on;
- plot(Am_freqs_kHz*1000, real(d), 'c')
- hold off;
-
- % decimated in time
-
- %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
-
- if plot_all_masks
- mask_sample_freqs_kHz = (1:L)*Wo*4/pi;
- for m=1:L
- maskdB = schroeder(m*Wo*4/pi, mask_sample_freqs_kHz) + AmdB(m);
- plot(mask_sample_freqs_kHz*1000, maskdB, "k--");
- end
- end
-
- hold off;
-
- if phase_stuff
-
- [phase Sdb s Aw] = determine_phase(model, f, ak(f,:));
- figure(3)
- subplot(211)
- plot(Sdb)
- title('Mag (dB)');
- subplot(212)
- plot(phase(1:256))
- hold on;
- plot(angle(Aw(1:256))+0.5,'g')
- hold off;
- title('Phase (rads)');
- figure(4)
- plot(s)
- end
+ figure(5)
+ clf
+ stem(dk_)
- % interactive menu
-
- printf("\rframe: %d menu: n-next b-back a-Am ", f);
- if freq_quant
- printf("F");
- else
- printf("f");
- end
- printf("-freq_quant ");
- if amp_quant
- printf("M");
- else
- printf("m");
- end
- printf("-amp_quant q-quit");
+ % interactive menu ------------------------------------------
+ printf("\rframe: %d menu: n-next b-back q-quit", f);
fflush(stdout);
k = kbhit();
+
if (k == 'n')
f = f + 1;
endif
if (k == 'b')
f = f - 1;
endif
- if k == 'm'
- if amp_quant == 0
- amp_quant = 3;
- else
- amp_quant = 0;
- end
- end
- if k == 'f'
- if freq_quant == 0
- freq_quant = 3;
- else
- freq_quant = 0;
- end
- end
- if k == 'a'
- if plot_spectrum == 0
- plot_spectrum = 1;
- else
- plot_spectrum = 0;
- end
- end
until (k == 'q')
printf("\n");
endfunction
+