1% tnewamp1.m
2%
3% Copyright David Rowe 2017
4% This program is distributed under the terms of the GNU General Public License
5% Version 2
6
7#{
8
9  Octave script to compare Octave and C versions of newamp1 processing, in order to test C port.
10
11  c2sim -> dump files -> $ ../build_linux/unittest/tnewamp1 -> octave:1> tnewamp1
12  Usage:
13
14    1/ build codec2 with -DDUMP - see codec2-dev/README
15
16    2/ Generate dump files using c2sim (just need to do this once)
17       $ cd codec2-dev/build_linux/src
18       $ ./c2sim ../../raw/hts1a.raw --phase0 --postfilter --dump hts1a --lpc 10 --dump_pitch_e hts1a_pitche.txt
19
20    3/ Run C version which generates a file of Octave test vectors as output:
21
22      $ cd codec2-dev/build_linux/unittest
23      $ ./tnewamp1 ../../raw/hts1a.raw
24
25    4/ Run Octave script to generate Octave test vectors and compare with C.
26
27      octave:1> tnewamp1("../build_linux/src/hts1a")
28
29    5/ Optionally listen to output
30
31     ~/codec2-dev/build_linux/src$ ./c2sim ../../raw/hts1a.raw --phase0 --postfilter \
32                                   --amread hts1a_am.out --hmread hts1a_hm.out \
33                                   --Woread hts1a_Wo.out --hand_voicing hts1a_v.txt -o - \
34                                      | play -q -t raw -r 8000 -s -2 -
35#}
36
37function tnewamp1(input_prefix, path_to_unittest="../build_linux/unittest/")
38  printf("starting tnewamp1.c input_prefix: %s\n", input_prefix);
39
40  visible_flag = 'off';
41  newamp_700c;
42  autotest;
43  more off;
44
45  max_amp = 80;
46  postfilter = 0;   % optional postfiler that runs on Am, not used atm
47  synth_phase = 1;
48
49  if nargin == 1
50    output_prefix = input_prefix;
51  end
52  model_name = strcat(input_prefix,"_model.txt");
53  model = load(model_name);
54  [frames nc] = size(model);
55
56  voicing_name = strcat(input_prefix,"_pitche.txt");
57  voicing = zeros(1,frames);
58
59  if exist(voicing_name, "file") == 2
60    pitche = load(voicing_name);
61    voicing = pitche(:, 3);
62  end
63
64  % Load in C vectors and compare -----------------------------------------
65
66  load(sprintf("%s/tnewamp1_out.txt", path_to_unittest));
67
68  K = 20;
69  [frames tmp] = size(rate_K_surface_c);
70  [rate_K_surface sample_freqs_kHz] = resample_const_rate_f_mel(model(1:frames,:), K);
71
72  melvq;
73  load train_120_1.txt; load train_120_2.txt;
74  train_120_vq(:,:,1)= train_120_1; train_120_vq(:,:,2)= train_120_2; m=5;
75  m=5;
76
77  eq = zeros(1,K);
78  for f=1:frames
79    mean_f(f) = mean(rate_K_surface(f,:));
80    rate_K_surface_no_mean(f,:) = rate_K_surface(f,:) - mean_f(f);
81    [rate_K_vec eq] = front_eq(rate_K_surface_no_mean(f,:), eq);
82    rate_K_surface_no_mean(f,:) = rate_K_vec;
83  end
84
85  [res rate_K_surface_no_mean_ ind] = mbest(train_120_vq, rate_K_surface_no_mean, m);
86
87  for f=1:frames
88    rate_K_surface_no_mean_(f,:) = post_filter(rate_K_surface_no_mean_(f,:), sample_freqs_kHz, 1.5);
89  end
90
91  rate_K_surface_ = zeros(frames, K);
92  interpolated_surface_ = zeros(frames, K);
93  energy_q = create_energy_q;
94  M = 4;
95  for f=1:frames
96    [mean_f_ indx] = quantise(energy_q, mean_f(f));
97    indexes(f,3) = indx - 1;
98    rate_K_surface_(f,:) = rate_K_surface_no_mean_(f,:) + mean_f_;
99  end
100
101  % simulated decoder
102  % break into segments of M frames.  We have 2 samples spaced M apart
103  % and interpolate the rest.
104
105  Nfft_phase = 128;  % note this needs to be 512 (FFT_ENC in codec2 if using --awread)
106                     % with --hmread 128 is preferred as less memory/CPU
107  model_ = zeros(frames, max_amp+2);
108  voicing_ = zeros(1,frames);
109  Aw = zeros(frames, Nfft_phase);
110  H = zeros(frames, max_amp);
111  model_(1,1) = Wo_left = 2*pi/100;
112  voicing_left = 0;
113  left_vec = zeros(1,K);
114
115  % decoder runs on every M-th frame, 25Hz frame rate, offset at
116  % start is to minimise processing delay (thanks Jeroen!)
117
118  for f=M:M:frames
119
120    if voicing(f)
121      index = encode_log_Wo(model(f,1), 6);
122      if index == 0
123        index = 1;
124      end
125      model_(f,1) = decode_log_Wo(index, 6);
126    else
127      model_(f,1) = 2*pi/100;
128    end
129
130    Wo_right = model_(f,1);
131    voicing_right = voicing(f);
132    [Wo_ avoicing_] = interp_Wo_v(Wo_left, Wo_right, voicing_left, voicing_right);
133
134    #{
135    for i=1:4
136      fprintf(stderr, "  Wo: %4.3f L: %d v: %d\n", Wo_(i), floor(pi/Wo_(i)), avoicing_(i));
137    end
138    fprintf(stderr,"  rate_K_vec: ");
139    for i=1:5
140      fprintf(stderr,"%5.3f  ", rate_K_surface_(f,i));
141    end
142    fprintf(stderr,"\n");
143    #}
144
145    if f > M
146      model_(f-M:f-1,1) = Wo_;
147      voicing_(f-M:f-1) = avoicing_;
148      model_(f-M:f-1,2) = floor(pi ./ model_(f-M:f-1,1)); % calculate L for each interpolated Wo
149    end
150
151    right_vec = rate_K_surface_(f,:);
152
153    if f > M
154      sample_points = [f-M f];
155      resample_points = f-M:f-1;
156      for k=1:K
157        interpolated_surface_(resample_points,k) = interp_linear(sample_points, [left_vec(k) right_vec(k)], resample_points);
158      end
159
160      for k=f-M:f-1
161        model_(k,:) = resample_rate_L(model_(k,:), interpolated_surface_(k,:), sample_freqs_kHz);
162        Aw(k,:) = determine_phase(model_, k, Nfft_phase);
163        for m=1:model_(k,2)
164          b = round(m*model_(k,1)*Nfft_phase/(2*pi));  % map harmonic centre to DFT bin
165          H(k,m) = exp(j*Aw(k, b+1));
166        end
167      end
168
169   end
170
171   % update for next time
172
173   Wo_left = Wo_right;
174   voicing_left = voicing_right;
175   left_vec = right_vec;
176
177  end
178
179  f = figure(1); clf;
180  mesh(angle(H));
181  f = figure(2); clf;
182  mesh(angle(H_c(:,1:max_amp)));
183  f = figure(3); clf;
184  mesh(abs(H - H_c(:,1:max_amp)));
185
186  passes = 0; tests = 0;
187  passes += check(eq, eq_c, 'Equaliser', 0.01); tests++;
188  passes += check(rate_K_surface, rate_K_surface_c, 'rate_K_surface', 0.01); tests++;
189  passes += check(mean_f, mean_c, 'mean', 0.01); tests++;
190  passes += check(rate_K_surface_, rate_K_surface__c, 'rate_K_surface_', 0.01); tests++;
191  passes += check(interpolated_surface_, interpolated_surface__c, 'interpolated_surface_', 0.01); tests++;
192  passes += check(model_(:,1), model__c(:,1), 'interpolated Wo_', 0.001);  tests++;
193  passes += check(voicing_, voicing__c, 'interpolated voicing'); tests++;
194  passes += check(model_(:,3:max_amp+2), model__c(:,3:max_amp+2), 'rate L Am surface ', 0.1); tests++;
195  passes += check(H, H_c(:,1:max_amp), 'phase surface'); tests++;
196  printf("passes: %d fails: %d\n", passes, tests - passes);
197
198  #{
199  % Save to disk to check synthesis is OK with c2sim
200
201  output_prefix = input_prefix;
202  Am_out_name = sprintf("%s_am.out", output_prefix);
203  fam  = fopen(Am_out_name,"wb");
204
205  Wo_out_name = sprintf("%s_Wo.out", output_prefix);
206  fWo  = fopen(Wo_out_name,"wb");
207
208  Aw_out_name = sprintf("%s_aw.out", output_prefix);
209  faw = fopen(Aw_out_name,"wb");
210
211  Hm_out_name = sprintf("%s_hm.out", output_prefix);
212  fhm = fopen(Hm_out_name,"wb");
213
214  printf("Generating files for c2sim: ");
215  for f=1:frames
216    printf(".", f);
217    Wo = model_(f,1);
218    L = min([model_(f,2) max_amp-1]);
219    Am = model_(f,3:(L+2));
220
221    Am_ = zeros(1,2*max_amp);
222    Am_(2:L) = Am(1:L-1);
223
224    fwrite(fam, Am_, "float32");
225    fwrite(fWo, Wo, "float32");
226
227    % Note we send opposite phase as c2sim expects phase of LPC
228    % analysis filter, just a convention based on historical
229    % development of Codec 2
230
231    Aw1 = zeros(1, Nfft_phase*2);
232    Aw1(1:2:Nfft_phase*2) = cos(Aw(f,:));
233    Aw1(2:2:Nfft_phase*2) = -sin(Aw(f,:));
234    fwrite(faw, Aw1, "float32");
235
236    Hm = zeros(1, 2*2*max_amp);
237    for m=1:L
238        Hm(2*m+1) = real(H(f,m));
239        Hm(2*m+2) = imag(H(f,m));
240    end
241    fwrite(fhm, Hm, "float32");
242  end
243
244  fclose(fam); fclose(fWo); fclose(faw); fclose(fhm);
245
246  v_out_name = sprintf("%s_v.txt", output_prefix);
247  fv  = fopen(v_out_name,"wt");
248  for f=1:length(voicing__c)
249    fprintf(fv,"%d\n", voicing__c(f));
250  end
251  fclose(fv);
252  #}
253
254endfunction
255
256
257