1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
|
#include "pv.h"
/* S is a spectrum in rfft format, i.e., it contains N real values
arranged as real followed by imaginary values, except for first
two values, which are real parts of 0 and Nyquist frequencies;
convert first changes these into N/2+1 PAIRS of magnitude and
phase values to be stored in output array C; the phases are then
unwrapped and successive phase differences are used to compute
estimates of the instantaneous frequencies for each phase vocoder
analysis channel; decimation rate D and sampling rate R are used
to render these frequency values directly in Hz. */
void convert(float *S, float *C, int N2, float *lastphase, float fundamental, float factor )
{
#if 1
float phase,phasediff;
int even,odd;
float a,b;
int i;
for ( i = 0; i <= N2; i++ ) {
odd = ( even = i<<1 ) + 1;
a = ( i == N2 ? S[1] : S[even] );
b = ( i == 0 || i == N2 ? 0. : S[odd] );
C[even] = hypot( a, b );
if ( C[even] == 0. )
phasediff = 0.;
else {
phase = -atan2( b, a );
phasediff = fmod(phase - lastphase[i] + (PV_2PI+PV_PI), PV_2PI)-PV_PI;
lastphase[i] = phase;
}
C[odd] = phasediff*factor + i*fundamental;
}
#else
float phase,
phasediff;
int real,
imag,
amp,
freq;
float a,
b;
int i;
float myTWOPI, myPI;
myTWOPI = 8.*atan(1.);
myPI = 4.*atan(1.);
for ( i = 0; i <= N2; i++ ) {
imag = freq = ( real = amp = i<<1 ) + 1;
a = ( i == N2 ? S[1] : S[real] );
b = ( i == 0 || i == N2 ? 0. : S[imag] );
C[amp] = hypot( a, b );
if ( C[amp] == 0. )
phasediff = 0.;
else {
phasediff = ( phase = -atan2( b, a ) ) - lastphase[i];
lastphase[i] = phase;
// TG: DANGEROUS!!!! (and slow, if lastphase not correctly initialized)
while ( phasediff > myPI )
phasediff -= myTWOPI;
while ( phasediff < -myPI )
phasediff += myTWOPI;
}
C[freq] = phasediff*factor + i*fundamental;
/*
if( i > 8 && i < 12 ) {
fprintf(stderr,"convert freq %d: %f\n",i, C[freq]);
}
*/
}
#endif
}
void unconvert(float *C, float *S, int N2, float *lastphase, float fundamental, float factor )
{
#if 1
int i,even,odd;
float mag,phase;
for ( i = 0; i <= N2; i++ ) {
odd = ( even = i<<1 ) + 1;
mag = C[even];
lastphase[i] += C[odd] - i*fundamental;
phase = lastphase[i]*factor;
if(i != N2) {
S[even] = mag*cos( phase );
S[odd] = -mag*sin( phase );
}
else
S[1] = mag*cos( phase );
}
#else
int i,
real,
imag,
amp,
freq;
float mag,
phase;
for ( i = 0; i <= N2; i++ ) {
imag = freq = ( real = amp = i<<1 ) + 1;
if ( i == N2 )
real = 1;
mag = C[amp];
lastphase[i] += C[freq] - i*fundamental;
phase = lastphase[i]*factor;
S[real] = mag*cos( phase );
if ( i != N2 )
S[imag] = -mag*sin( phase );
/*
if( i == 10 ) {
fprintf(stderr,"unconvert: amp: %f freq: %f funda %f fac %f\n", C[amp],C[freq],fundamental,factor);
}
*/
}
#endif
}
|
