1% vq_700c.m 2% David Rowe May 2019 3% 4% Researching Codec 2 700C VQ equaliser ideas 5% See also scripts/train_700c_quant.sh, tnewamp1.m 6 7melvq; 8newamp_700c; 9 10% general purpose plot function for looking at averages of K-band 11% sequences in scripts dir and VQs: 12% vq_700c_plots({"hts2a.f32" "vk5qi.f32" "train_120_1.txt"}) 13 14function vq_700c_plots(fn_array) 15 K = 20; rate_K_sample_freqs_kHz = mel_sample_freqs_kHz(K); 16 freq_Hz = rate_K_sample_freqs_kHz * 1000; 17 18 figure(1); clf; hold on; axis([200 4000 40 90]); title('Max Hold'); 19 figure(2); clf; hold on; axis([200 4000 0 40]); title('Average'); 20 21 for i=1:length(fn_array) 22 [dir name ext] = fileparts(fn_array{i}); 23 if strcmp(ext, ".f32") 24 % f32 feature file 25 fn = sprintf("../build_linux/%s%s", name, ext) 26 bands = load_f32(fn , K); 27 else 28 % text file (e.g. existing VQ) 29 bands = load(fn_array{i}); 30 end 31 % for max hold: break into segments of Nsec, find max, average maximums 32 % this avoids very rare global peaks setting the max 33 Nsec = 10; Tframe = 0.01; frames_per_seg = Nsec/Tframe 34 Nsegs = floor(length(bands)/frames_per_seg) 35 max_holds = zeros(Nsegs, K); 36 if Nsegs == 0 37 max_holds = max(bands) 38 else 39 for s=1:Nsegs 40 st = (s-1)*frames_per_seg+1; en = st + frames_per_seg - 1; 41 max_holds(s,:) = max(bands(st:en,:)); 42 end 43 max_holds = mean(max_holds); 44 end 45 figure(1); plot(freq_Hz, max_holds, '+-', 'linewidth', 2); 46 figure(2); plot(freq_Hz, mean(bands), '+-', 'linewidth', 2); 47 end 48 figure(1); legend(fn_array); grid; xlabel('Freq (Hz)'); ylabel('Amp dB'); 49 figure(2); legend(fn_array); grid; xlabel('Freq (Hz)'); ylabel('Amp dB'); 50endfunction 51 52 53% limit mean of each vector to between lower_lim and upper_lim 54function vout = limit_vec(vin, lower_lim, upper_lim) 55 m = mean(vin'); 56 vout = zeros(size(vin)); 57 for i=1:length(vin) 58 vec_no_mean = vin(i,:) - m(i); 59 if m(i) < lower_lim 60 m(i) = lower_lim; 61 end 62 if m(i) > upper_lim 63 m(i) = upper_lim; 64 end 65 vout(i,:) = vec_no_mean + m(i); 66 end 67endfunction 68 69 70% single stage vq a target matrix 71function errors = vq_targets(vq, targets) 72 errors = []; 73 for i=1:length(targets) 74 [mse_list index_list] = search_vq(vq, targets(i,:), 1); 75 error = targets(i,:) - vq(index_list(1),:); 76 errors = [errors; error]; 77 end 78endfunction 79 80 81% single stage vq a target matrix with adaptive EQ, this didn't work 82 83function [errors eqs] = vq_targets_adap_eq(vq, targets, eqs) 84 errors = []; gain=0.02; 85 eq = eqs(end,:); 86 for i=1:length(targets) 87 t = targets(i,:) - eq; 88 mean(t) 89 %t -= mean(t); 90 [mse_list index_list] = search_vq(vq, t, 1); 91 error = t - vq(index_list(1),:); 92 eq = (1-gain)*eq + gain*error; 93 errors = [errors; error]; eqs = [eqs; eq]; 94 end 95endfunction 96 97 98% single stage vq a target matrix with block adaptive EQ, this works 99% well with nblock == 10 100 101function [errors eq] = vq_targets_block_eq(vq, targets, eq, nblock) 102 errors = []; n = 0; [tmp K] = size(vq); error_eq = zeros(1,K); gain=0.20; 103 for i=1:length(targets) 104 t = targets(i,:) - eq; 105 [mse_list index_list] = search_vq(vq, t, 1); 106 error = t - vq(index_list(1),:); 107 error_eq += error; 108 errors = [errors; error]; 109 n++; 110 if n == nblock 111 eq = 0.99*eq + gain*error_eq/nblock; 112 n = 0; error_eq = zeros(1,K); 113 end 114 end 115endfunction 116 117 118% two stage mbest VQ a target matrix 119 120function [errors targets_] = vq_targets2(vq1, vq2, targets) 121 vqset(:,:,1)= vq1; vqset(:,:,2)=vq2; m=5; 122 [errors targets_] = mbest(vqset, targets, m); 123endfunction 124 125 126% two stage mbest VQ a target matrix, with adap_eq 127 128function [errors targets_ eq] = vq_targets2_adap_eq(vq1, vq2, targets, eq) 129 vqset(:,:,1)= vq1; vqset(:,:,2)=vq2; m=5; gain=0.02; 130 errors = []; targets_ = []; 131 for i=1:length(targets) 132 t = targets(i,:)-eq; 133 t -= mean(t')'; 134 [error target_ indexes] = mbest(vqset, t, m); 135 % use first stage VQ as error driving adaptive EQ 136 eq_error = t - vq1(indexes(1),:); 137 eq = (1-gain)*eq + gain*eq_error; 138 errors = [errors; error]; targets_ = [targets_; target_]; 139 end 140endfunction 141 142 143% Given target and vq matrices, estimate eq via two metrics. First 144% metric seems to work best. Both uses first stage VQ error for EQ 145 146function [eq1 eq2] = est_eq(vq, targets) 147 [ntargets K] = size(targets); 148 [nvq K] = size(vq); 149 150 eq1 = zeros(1,K); eq2 = zeros(1,K); 151 for i=1:length(targets) 152 [mse_list index_list] = search_vq(vq, targets(i,:), 1); 153 154 % eq metric 1: average of error for best VQ entry 155 eq1 += targets(i,:) - vq(index_list(1),:); 156 157 % eq metric 2: average of error across all VQ entries 158 for j=1:nvq 159 eq2 += targets(i,:) - vq(j,:); 160 end 161 end 162 163 eq1 /= ntargets; 164 eq2 /= (ntargets*nvq); 165endfunction 166 167function [targets e] = load_targets(fn_target_f32) 168 nb_features = 41; 169 K = 20; 170 171 % .f32 files are in scripts directory, first K values rate_K_no_mean vectors 172 [dir name ext] = fileparts(fn_target_f32); 173 fn = sprintf("../script/%s_feat.f32", name); 174 feat = load_f32(fn, nb_features); 175 e = feat(:,1); 176 targets = feat(:,2:K+1); 177endfunction 178 179% rather simple EQ in front of VQ 180 181function [eqs ideal] = est_eq_front(targets) 182 [tmp K] = size(targets); 183 ideal = [ 8 10 12 14 14*ones(1,K-1-4) -20]; 184 eq = zeros(1,K); gain = 0.02; 185 eqs = []; 186 for i=1:length(targets) 187 update = targets(i,:) - ideal; 188 eq = (1-gain)*eq + gain*update; 189 eq(find(eq < 0)) = 0; 190 eqs = [eqs; eq]; 191 end 192endfunction 193 194function table_across_samples 195 K = 20; 196 197 % VQ is in .txt file in this directory, we have two to choose from. train_120 is the Codec 2 700C VQ, 198 % train_all_speech was trained up from a different, longer database, as a later exercise 199 vq_name = "train_120"; 200 #vq_name = "train_all_speech"; 201 vq1 = load(sprintf("%s_1.txt", vq_name)); 202 vq2 = load(sprintf("%s_2.txt", vq_name)); 203 204 printf("----------------------------------------------------------------------------------\n"); 205 printf("Sample Initial vq1 vq1_eq2 vq1_eq2 vq2 vq2_eq1 vq2_eq2 \n"); 206 printf("----------------------------------------------------------------------------------\n"); 207 208 fn_targets = { "cq_freedv_8k_lfboost" "cq_freedv_8k_hfcut" "cq_freedv_8k" "hts1a" "hts2a" "cq_ref" "ve9qrp_10s" "vk5qi" "c01_01_8k" "ma01_01"}; 209 #fn_targets = {"cq_freedv_8k_lfboost"}; 210 figs=1; 211 for i=1:length(fn_targets) 212 213 % load target and estimate eq 214 [targets e] = load_targets(fn_targets{i}); 215 eq1 = est_eq(vq1, targets); 216 eq2s = est_eq_front(targets); 217 % for these simulation uses fixed EQ sample, rather than letting it vary frame by frame 218 eq2 = eq2s(end,:); 219 220 % first stage VQ ----------------- 221 222 errors1 = vq_targets(vq1, targets); 223 errors1_eq1 = vq_targets(vq1, targets-eq1); 224 errors1_eq2 = vq_targets(vq1, targets-eq2); 225 226 % two stage mbest VQ -------------- 227 228 [errors2 targets_] = vq_targets2(vq1, vq2, targets); 229 [errors2_eq1 targets_eq1_] = vq_targets2(vq1, vq2, targets-eq1); 230 [errors2_eq2 targets_eq2_] = vq_targets2(vq1, vq2, targets-eq2); 231 232 % save to .f32 files for listening tests 233 if strcmp(vq_name,"train_120") 234 save_f32(sprintf("../script/%s_vq2.f32", fn_targets{i}), targets_); 235 save_f32(sprintf("../script/%s_vq2_eq1.f32", fn_targets{i}), targets_eq1_); 236 save_f32(sprintf("../script/%s_vq2_eq2.f32", fn_targets{i}), targets_eq2_); 237 else 238 save_f32(sprintf("../script/%s_vq2_as.f32", fn_targets{i}), targets_); 239 save_f32(sprintf("../script/%s_vq2_as_eq.f32", fn_targets{i}), targets_eq_); 240 end 241 printf("%-21s %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\n", fn_targets{i}, 242 var(targets(:)), var(errors1(:)), var(errors1_eq1(:)), var(errors1_eq2(:)), 243 var(errors2(:)), var(errors2_eq1(:)), var(errors2_eq2(:))); 244 245 figure(figs++); clf; 246 %plot(var(errors2'),'b;vq2;'); hold on; plot(var(errors2_eq1'),'g;vq2_eq1;'); plot(var(errors2_eq2'),'r;vq2_eq2;'); hold off; 247 plot(mean(targets),'g;mean(targets);'); hold on; plot(mean(vq1),'g;mean(vq1);'); plot(eq2,'r;eq2;'); hold off; 248 title(fn_targets{i}); axis([1 K -20 30]); 249 end 250endfunction 251 252 253% interactve, menu driven frame by frame plots 254 255function interactive(fn_vq_txt, fn_target_f32) 256 K = 20; 257 vq = load("train_120_1.txt"); 258 [targets e] = load_targets(fn_target_f32); 259 eq1 = est_eq(vq, targets); 260 261 [errors1_eq2 eqs2] = vq_targets_adap_eq(vq, targets, zeros(1,K)); 262 [errors1_eq2 eqs2] = vq_targets_adap_eq(vq, targets, eqs2(end,:)); 263 eq2 = eqs2(end,:); 264 265 figure(1); clf; 266 mesh(e+targets) 267 figure(2); clf; 268 plot(eq1,'b;eq1;') 269 hold on; 270 plot(mean(targets),'c;mean(targets);'); plot(eq2,'g;eq2;'); 271 hold off; 272 figure(3); clf; mesh(eqs2); title('eq2 evolving') 273 274 % enter single step loop 275 f = 20; neq = 0; eq=zeros(1,K); 276 do 277 figure(4); clf; 278 t = targets(f,:) - eq; 279 [mse_list index_list] = search_vq(vq, t, 1); 280 error = t - vq(index_list(1),:); 281 plot(e(f)+t,'b;target;'); 282 hold on; 283 plot(e(f)+vq(index_list,:),'g;vq;'); 284 plot(error,'r;error;'); 285 plot(eq,'c;eq;'); 286 plot([1 K],[e(f) e(f)],'--') 287 hold off; 288 axis([1 K -20 80]) 289 % interactive menu 290 291 printf("\r f: %2d eq: %d ind: %3d var: %3.1f menu: n-next b-back e-eq q-quit", f, neq, index_list(1), var(error)); 292 fflush(stdout); 293 k = kbhit(); 294 295 if k == 'n' f+=1; end 296 if k == 'e' 297 neq++; 298 end 299 if neq == 3 neq = 0; end 300 if neq == 0 eq = zeros(1,K); end 301 if neq == 1 eq = eq1; end 302 if neq == 2 eq = eqs2(f,:); end 303 if k == 'b' f-=1; end 304 until (k == 'q') 305 printf("\n"); 306endfunction 307 308 309% Experiment to test iterative approach of block update and remove 310% mean (ie frame energy), shows some promise at reducing HF energy 311% over several iterations while not affecting already good samples 312 313function experiment_iterate_block(fn_vq_txt, fn_target_f32) 314 K = 20; 315 vq = load("train_120_1.txt"); 316 [targets e] = load_targets(fn_target_f32); 317 318 figure(1); clf; 319 plot(mean(targets),'b;mean(targets);'); 320 hold on; 321 plot(mean(vq), 'g;mean(vq);'); 322 figure(2); clf; hold on; 323 eq = zeros(1,K); 324 for i=1:3 325 [errors eq] = vq_targets_block_eq(vq, targets, eq, 10); 326 figure(1); plot(mean(targets-eq)); 327 figure(2); plot(eq); 328 printf("i: %d %6.2f\n", i, var(errors(:))) 329 end 330endfunction 331 332% Experiment to test EQ of input (before) VQ. We set a threshold on 333% when to equalise, so we don't upset already flat-ish samples. This 334% is the algorithm used for C at the time of writing (newamp1.c, newamp_700c.m) 335 336function experiment_front_eq(fn_vq_txt, fn_target_f32) 337 K = 20; 338 vq = load("train_120_1.txt"); 339 [targets e] = load_targets(fn_target_f32); 340 341 [eqs ideal] = est_eq_front(targets); 342 343 figure(1); clf; 344 plot(mean(targets),'b;mean(targets);'); 345 hold on; 346 plot(ideal, 'g;ideal;'); 347 plot(eqs(end,:), 'r;eq;'); 348 plot(mean(targets)-eqs(end,:), 'c;equalised;'); 349 plot(mean(vq),'b--;mean(vq);'); 350 hold off; 351 figure(2); clf; mesh(eqs(1:100,:)); title('EQ weights over time'); 352 ylabel('Time (frames'); xlabel('Freq (mel)'); 353endfunction 354 355more off 356 357% choose one of these to run first 358% You'll need to run scripts/train_700C_quant.sh first to generate the .f32 files 359 360%interactive("train_120_1.txt", "cq_freedv_8k_lfboost.f32") 361%table_across_samples; 362%vq_700c_plots({"all_speech_8k.f32" "all_speech_8k_hp300.f32" "dev-clean-8k.f32" "train_8k.f32" } ) 363%vq_700c_plots({"ve9qrp_10s.f32" "cq_freedv_8k_lfboost.f32" "cq_ref.f32" "hts1a.f32" "vk5qi.f32"}) 364%experiment_iterate_block("train_120_1.txt", "ve9qrp_10s.f32") 365%experiment_iterate_block("train_120_1.txt", "cq_freedv_8k_lfboost.f32") 366%experiment_front_eq("train_120_1.txt", "cq_freedv_8k_lfboost.f32") 367