summaryrefslogtreecommitdiff
path: root/phase.c
diff options
context:
space:
mode:
Diffstat (limited to 'phase.c')
-rw-r--r--phase.c242
1 files changed, 114 insertions, 128 deletions
diff --git a/phase.c b/phase.c
index e486613..dec8793 100644
--- a/phase.c
+++ b/phase.c
@@ -25,18 +25,19 @@
25 along with this program; if not,see <http://www.gnu.org/licenses/>. 25 along with this program; if not,see <http://www.gnu.org/licenses/>.
26*/ 26*/
27 27
28#include "defines.h"
29#include "phase.h" 28#include "phase.h"
30#include "kiss_fft.h"
31#include "comp.h"
32#include "comp_prim.h"
33#include "sine.h"
34 29
35#include <assert.h> 30#include <assert.h>
36#include <ctype.h> 31#include <ctype.h>
37#include <math.h> 32#include <math.h>
38#include <string.h>
39#include <stdlib.h> 33#include <stdlib.h>
34#include <string.h>
35
36#include "comp.h"
37#include "comp_prim.h"
38#include "defines.h"
39#include "kiss_fft.h"
40#include "sine.h"
40 41
41/*---------------------------------------------------------------------------*\ 42/*---------------------------------------------------------------------------*\
42 43
@@ -47,25 +48,23 @@
47 48
48\*---------------------------------------------------------------------------*/ 49\*---------------------------------------------------------------------------*/
49 50
50void sample_phase(MODEL *model, 51void sample_phase(MODEL *model, COMP H[],
51 COMP H[], 52 COMP A[] /* LPC analysis filter in freq domain */
52 COMP A[] /* LPC analysis filter in freq domain */ 53) {
53 ) 54 int m, b;
54{ 55 float r;
55 int m, b;
56 float r;
57 56
58 r = TWO_PI/(FFT_ENC); 57 r = TWO_PI / (FFT_ENC);
59 58
60 /* Sample phase at harmonics */ 59 /* Sample phase at harmonics */
61 60
62 for(m=1; m<=model->L; m++) { 61 for (m = 1; m <= model->L; m++) {
63 b = (int)(m*model->Wo/r + 0.5); 62 b = (int)(m * model->Wo / r + 0.5);
64 H[m] = cconj(A[b]); /* synth filter 1/A is opposite phase to analysis filter */ 63 H[m] =
65 } 64 cconj(A[b]); /* synth filter 1/A is opposite phase to analysis filter */
65 }
66} 66}
67 67
68
69/*---------------------------------------------------------------------------*\ 68/*---------------------------------------------------------------------------*\
70 69
71 phase_synth_zero_order() 70 phase_synth_zero_order()
@@ -158,64 +157,56 @@ void sample_phase(MODEL *model,
158\*---------------------------------------------------------------------------*/ 157\*---------------------------------------------------------------------------*/
159 158
160void phase_synth_zero_order( 159void phase_synth_zero_order(
161 int n_samp, 160 int n_samp, MODEL *model,
162 MODEL *model, 161 float *ex_phase, /* excitation phase of fundamental */
163 float *ex_phase, /* excitation phase of fundamental */ 162 COMP H[] /* L synthesis filter freq domain samples */
164 COMP H[] /* L synthesis filter freq domain samples */ 163
165 164) {
166) 165 int m;
167{ 166 float new_phi;
168 int m; 167 COMP Ex[MAX_AMP + 1]; /* excitation samples */
169 float new_phi; 168 COMP A_[MAX_AMP + 1]; /* synthesised harmonic samples */
170 COMP Ex[MAX_AMP+1]; /* excitation samples */ 169
171 COMP A_[MAX_AMP+1]; /* synthesised harmonic samples */ 170 /*
172 171 Update excitation fundamental phase track, this sets the position
173 /* 172 of each pitch pulse during voiced speech. After much experiment
174 Update excitation fundamental phase track, this sets the position 173 I found that using just this frame's Wo improved quality for UV
175 of each pitch pulse during voiced speech. After much experiment 174 sounds compared to interpolating two frames Wo like this:
176 I found that using just this frame's Wo improved quality for UV 175
177 sounds compared to interpolating two frames Wo like this: 176 ex_phase[0] += (*prev_Wo+model->Wo)*N_SAMP/2;
178 177 */
179 ex_phase[0] += (*prev_Wo+model->Wo)*N_SAMP/2; 178
180 */ 179 ex_phase[0] += (model->Wo) * n_samp;
181 180 ex_phase[0] -= TWO_PI * floorf(ex_phase[0] / TWO_PI + 0.5);
182 ex_phase[0] += (model->Wo)*n_samp; 181
183 ex_phase[0] -= TWO_PI*floorf(ex_phase[0]/TWO_PI + 0.5); 182 for (m = 1; m <= model->L; m++) {
184 183 /* generate excitation */
185 for(m=1; m<=model->L; m++) { 184
186 185 if (model->voiced) {
187 /* generate excitation */ 186 Ex[m].real = cosf(ex_phase[0] * m);
188 187 Ex[m].imag = sinf(ex_phase[0] * m);
189 if (model->voiced) { 188 } else {
190 189 /* When a few samples were tested I found that LPC filter
191 Ex[m].real = cosf(ex_phase[0]*m); 190 phase is not needed in the unvoiced case, but no harm in
192 Ex[m].imag = sinf(ex_phase[0]*m); 191 keeping it.
193 } 192 */
194 else { 193 float phi = TWO_PI * (float)codec2_rand() / CODEC2_RAND_MAX;
195 194 Ex[m].real = cosf(phi);
196 /* When a few samples were tested I found that LPC filter 195 Ex[m].imag = sinf(phi);
197 phase is not needed in the unvoiced case, but no harm in 196 }
198 keeping it.
199 */
200 float phi = TWO_PI*(float)codec2_rand()/CODEC2_RAND_MAX;
201 Ex[m].real = cosf(phi);
202 Ex[m].imag = sinf(phi);
203 }
204
205 /* filter using LPC filter */
206 197
207 A_[m].real = H[m].real*Ex[m].real - H[m].imag*Ex[m].imag; 198 /* filter using LPC filter */
208 A_[m].imag = H[m].imag*Ex[m].real + H[m].real*Ex[m].imag;
209 199
210 /* modify sinusoidal phase */ 200 A_[m].real = H[m].real * Ex[m].real - H[m].imag * Ex[m].imag;
201 A_[m].imag = H[m].imag * Ex[m].real + H[m].real * Ex[m].imag;
211 202
212 new_phi = atan2f(A_[m].imag, A_[m].real+1E-12); 203 /* modify sinusoidal phase */
213 model->phi[m] = new_phi;
214 }
215 204
205 new_phi = atan2f(A_[m].imag, A_[m].real + 1E-12);
206 model->phi[m] = new_phi;
207 }
216} 208}
217 209
218
219/*---------------------------------------------------------------------------*\ 210/*---------------------------------------------------------------------------*\
220 211
221 FUNCTION....: mag_to_phase 212 FUNCTION....: mag_to_phase
@@ -230,60 +221,55 @@ void phase_synth_zero_order(
230 221
231\*---------------------------------------------------------------------------*/ 222\*---------------------------------------------------------------------------*/
232 223
233void mag_to_phase(float phase[], /* Nfft/2+1 output phase samples in radians */ 224void mag_to_phase(
234 float Gdbfk[], /* Nfft/2+1 postive freq amplitudes samples in dB */ 225 float phase[], /* Nfft/2+1 output phase samples in radians */
235 int Nfft, 226 float Gdbfk[], /* Nfft/2+1 positive freq amplitudes samples in dB */
236 codec2_fft_cfg fft_fwd_cfg, 227 int Nfft, codec2_fft_cfg fft_fwd_cfg, codec2_fft_cfg fft_inv_cfg) {
237 codec2_fft_cfg fft_inv_cfg 228 COMP Sdb[Nfft], c[Nfft], cf[Nfft], Cf[Nfft];
238 ) 229 int Ns = Nfft / 2 + 1;
239{ 230 int i;
240 COMP Sdb[Nfft], c[Nfft], cf[Nfft], Cf[Nfft]; 231
241 int Ns = Nfft/2+1; 232 /* install negative frequency components, 1/Nfft takes into
242 int i; 233 account kiss fft lack of scaling on ifft */
243 234
244 /* install negative frequency components, 1/Nfft takes into 235 Sdb[0].real = Gdbfk[0];
245 account kiss fft lack of scaling on ifft */ 236 Sdb[0].imag = 0.0;
246 237 for (i = 1; i < Ns; i++) {
247 Sdb[0].real = Gdbfk[0]; 238 Sdb[i].real = Sdb[Nfft - i].real = Gdbfk[i];
248 Sdb[0].imag = 0.0; 239 Sdb[i].imag = Sdb[Nfft - i].imag = 0.0;
249 for(i=1; i<Ns; i++) { 240 }
250 Sdb[i].real = Sdb[Nfft-i].real = Gdbfk[i]; 241
251 Sdb[i].imag = Sdb[Nfft-i].imag = 0.0; 242 /* compute real cepstrum from log magnitude spectrum */
252 } 243
253 244 codec2_fft(fft_inv_cfg, Sdb, c);
254 /* compute real cepstrum from log magnitude spectrum */ 245 for (i = 0; i < Nfft; i++) {
255 246 c[i].real /= (float)Nfft;
256 codec2_fft(fft_inv_cfg, Sdb, c); 247 c[i].imag /= (float)Nfft;
257 for(i=0; i<Nfft; i++) { 248 }
258 c[i].real /= (float)Nfft; 249
259 c[i].imag /= (float)Nfft; 250 /* Fold cepstrum to reflect non-min-phase zeros inside unit circle */
260 } 251
261 252 cf[0] = c[0];
262 /* Fold cepstrum to reflect non-min-phase zeros inside unit circle */ 253 for (i = 1; i < Ns - 1; i++) {
263 254 cf[i] = cadd(c[i], c[Nfft - i]);
264 cf[0] = c[0]; 255 }
265 for(i=1; i<Ns-1; i++) { 256 cf[Ns - 1] = c[Ns - 1];
266 cf[i] = cadd(c[i],c[Nfft-i]); 257 for (i = Ns; i < Nfft; i++) {
267 } 258 cf[i].real = 0.0;
268 cf[Ns-1] = c[Ns-1]; 259 cf[i].imag = 0.0;
269 for(i=Ns; i<Nfft; i++) { 260 }
270 cf[i].real = 0.0; 261
271 cf[i].imag = 0.0; 262 /* Cf = dB_magnitude + j * minimum_phase */
272 } 263
273 264 codec2_fft(fft_fwd_cfg, cf, Cf);
274 /* Cf = dB_magnitude + j * minimum_phase */ 265
275 266 /* The maths says we are meant to be using log(x), not 20*log10(x),
276 codec2_fft(fft_fwd_cfg, cf, Cf); 267 so we need to scale the phase to account for this:
277 268 log(x) = 20*log10(x)/scale */
278 /* The maths says we are meant to be using log(x), not 20*log10(x), 269
279 so we need to scale the phase to account for this: 270 float scale = (20.0 / logf(10.0));
280 log(x) = 20*log10(x)/scale */ 271
281 272 for (i = 0; i < Ns; i++) {
282 float scale = (20.0/logf(10.0)); 273 phase[i] = Cf[i].imag / scale;
283 274 }
284 for(i=0; i<Ns; i++) {
285 phase[i] = Cf[i].imag/scale;
286 }
287
288
289} 275}