pitch.c 5.91 KB
Newer Older
1
2
/* (C) 2007-2008 Jean-Marc Valin, CSIRO
*/
3
4
5
/**
   @file pitch.c
   @brief Pitch analysis
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
/*
   Redistribution and use in source and binary forms, with or without
   modification, are permitted provided that the following conditions
   are met:
   
   - Redistributions of source code must retain the above copyright
   notice, this list of conditions and the following disclaimer.
   
   - Redistributions in binary form must reproduce the above copyright
   notice, this list of conditions and the following disclaimer in the
   documentation and/or other materials provided with the distribution.
   
   - Neither the name of the Xiph.org Foundation nor the names of its
   contributors may be used to endorse or promote products derived from
   this software without specific prior written permission.
   
   THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
   ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
   LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
   A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR
   CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
   EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
   PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
   PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
   LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
   NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
   SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35
36
*/

37

38
39
40
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
41

42
43
44
45
/*#include "_kiss_fft_guts.h"
#include "kiss_fftr.h"*/
#include "kfft_single.h"

46
#include "pitch.h"
47
#include "psy.h"
48
#include "os_support.h"
49
#include "mathops.h"
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

kiss_fftr_cfg pitch_state_alloc(int max_lag)
{
   return kiss_fftr_alloc_celt_single(max_lag, 0, 0);
}

void pitch_state_free(kiss_fftr_cfg st)
{
   kiss_fft_free(st);
}

#ifdef FIXED_POINT
static void normalise16(celt_word16_t *x, int len, celt_word16_t val)
{
   int i;
   celt_word16_t maxval = 0;
   for (i=0;i<len;i++)
      maxval = MAX16(maxval, ABS16(x[i]));
   if (maxval > val)
   {
      int shift = 0;
      while (maxval > val)
      {
         maxval >>= 1;
         shift++;
      }
      if (shift==0)
         return;
      for (i=0;i<len;i++)
         x[i] = SHR16(x[i], shift);
   } else {
      int shift=0;
      if (maxval == 0)
         return;
      val >>= 1;
      while (maxval < val)
      {
         val >>= 1;
         shift++;
      }
      if (shift==0)
         return;
      for (i=0;i<len;i++)
         x[i] = SHL16(x[i], shift);
   }
}
#else
#define normalise16(x,len,val)
#endif

#define INPUT_SHIFT 15
101

102
void find_spectral_pitch(kiss_fftr_cfg fft, const struct PsyDecay *decay, const celt_sig_t *x, const celt_sig_t *y, const celt_word16_t *window, int overlap, int lag, int len, int C, int *pitch)
103
{
104
   int c, i;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
105
   celt_word32_t max_corr;
106
107
108
   VARDECL(celt_word16_t, X);
   VARDECL(celt_word16_t, Y);
   VARDECL(celt_mask_t, curve);
109
   int n2;
110
   int L2;
111
   const int *bitrev;
112
   SAVE_STACK;
113
   n2 = lag/2;
114
   L2 = len/2;
115
   ALLOC(X, lag, celt_word16_t);
116
   ALLOC(curve, n2, celt_mask_t);
117

118
   bitrev = fft->substate->bitrev;
119
   for (i=0;i<lag;i++)
120
      X[i] = 0;
121
   /* Sum all channels of the current frame and copy into X in bit-reverse order */
122
   for (c=0;c<C;c++)
123
   {
124
      for (i=0;i<L2;i++)
125
      {
126
127
         X[2*bitrev[i]] += SHR32(x[C*(2*i)+c],INPUT_SHIFT);
         X[2*bitrev[i]+1] += SHR32(x[C*(2*i+1)+c],INPUT_SHIFT);
128
129
      }
   }
130
131
   /* Applying the window in the bit-reverse domain. It's a bit weird, but it
      can help save memory */
132
133
   for (i=0;i<overlap/2;i++)
   {
134
135
136
137
      X[2*bitrev[i]] = MULT16_16_Q15(window[2*i], X[2*bitrev[i]]);
      X[2*bitrev[i]+1] = MULT16_16_Q15(window[2*i+1], X[2*bitrev[i]+1]);
      X[2*bitrev[L2-i-1]] = MULT16_16_Q15(window[2*i+1], X[2*bitrev[L2-i-1]]);
      X[2*bitrev[L2-i-1]+1] = MULT16_16_Q15(window[2*i], X[2*bitrev[L2-i-1]+1]);
138
   }
139
140
   normalise16(X, lag, 8192);
   /*for (i=0;i<lag;i++) printf ("%d ", X[i]);printf ("\n");*/
141
   /* Forward real FFT (in-place) */
142
143
144
   kf_work((kiss_fft_cpx*)X, NULL, 1,1, fft->substate->factors,fft->substate, 1, 1, 1);
   kiss_fftr_twiddles(fft,X);

145
   compute_masking(decay, X, curve, lag);
146
147

   /* Deferred allocation to reduce peak stack usage */
148
   ALLOC(Y, lag, celt_word16_t);
149
   for (i=0;i<lag;i++)
150
      Y[i] = 0;
151
   /* Sum all channels of the past audio and copy into Y in bit-reverse order */
152
   for (c=0;c<C;c++)
153
   {
154
      for (i=0;i<n2;i++)
155
      {
156
157
         Y[2*bitrev[i]] += SHR32(y[C*(2*i)+c],INPUT_SHIFT);
         Y[2*bitrev[i]+1] += SHR32(y[C*(2*i+1)+c],INPUT_SHIFT);
158
159
      }
   }
160
   normalise16(Y, lag, 8192);
161
   /* Forward real FFT (in-place) */
162
163
164
   kf_work((kiss_fft_cpx*)Y, NULL, 1,1, fft->substate->factors,fft->substate, 1, 1, 1);
   kiss_fftr_twiddles(fft,Y);

165
   /* Compute cross-spectrum using the inverse masking curve as weighting */
166
   for (i=1;i<n2;i++)
167
   {
168
      celt_word16_t n;
Jean-Marc Valin's avatar
Jean-Marc Valin committed
169
      celt_word32_t tmp;
170
      /*printf ("%d %d ", X[2*i]*X[2*i]+X[2*i+1]*X[2*i+1], Y[2*i]*Y[2*i]+Y[2*i+1]*Y[2*i+1]);*/
171
172
      /*n = DIV32_16(Q15ONE,celt_sqrt(EPSILON+curve[i]));*/
      n = ROUND(celt_rcp(celt_sqrt(EPSILON+curve[i])),16);
173
      /*printf ("%f ", n);*/
174
      tmp = X[2*i];
175
176
      X[2*i] = MULT16_32_Q15(n, ADD32(MULT16_16(X[2*i  ],Y[2*i  ]), MULT16_16(X[2*i+1],Y[2*i+1])));
      X[2*i+1] = MULT16_32_Q15(n, SUB32(MULT16_16(tmp,Y[2*i+1]), MULT16_16(X[2*i+1],Y[2*i  ])));
177
   }
178
   /*printf ("\n");*/
179
   X[0] = X[1] = 0;
180
181
   /*for (i=0;i<lag;i++) printf ("%d ", X[i]);printf ("\n");*/
   normalise16(X, lag, 50);
182
   /* Inverse half-complex to real FFT gives us the correlation */
183
   kiss_fftri(fft, X, Y);
184
   
185
   /* The peak in the correlation gives us the pitch */
Jean-Marc Valin's avatar
Jean-Marc Valin committed
186
   max_corr=-VERY_LARGE32;
187
   *pitch = 0;
188
189
   for (i=0;i<lag-len;i++)
   {
190
      /*printf ("%f ", xx[i]);*/
191
      if (Y[i] > max_corr)
192
193
      {
         *pitch = i;
194
         max_corr = Y[i];
195
196
      }
   }
197
   RESTORE_STACK;
198
}