qm-dsp  1.8
FFT.cpp
Go to the documentation of this file.
1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
2 
3 /*
4  QM DSP Library
5 
6  Centre for Digital Music, Queen Mary, University of London.
7  This file is based on Don Cross's public domain FFT implementation.
8 */
9 
10 #include "FFT.h"
11 
12 #include "maths/MathUtilities.h"
13 
14 #include "kiss_fft.h"
15 #include "kiss_fftr.h"
16 
17 #include <cmath>
18 
19 #include <iostream>
20 
21 #include <stdexcept>
22 
23 class FFT::D
24 {
25 public:
26  D(int n) : m_n(n) {
27  m_planf = kiss_fft_alloc(m_n, 0, NULL, NULL);
28  m_plani = kiss_fft_alloc(m_n, 1, NULL, NULL);
29  m_kin = new kiss_fft_cpx[m_n];
30  m_kout = new kiss_fft_cpx[m_n];
31  }
32 
33  ~D() {
34  kiss_fft_free(m_planf);
35  kiss_fft_free(m_plani);
36  delete[] m_kin;
37  delete[] m_kout;
38  }
39 
40  void process(bool inverse,
41  const double *ri,
42  const double *ii,
43  double *ro,
44  double *io) {
45 
46  for (int i = 0; i < m_n; ++i) {
47  m_kin[i].r = ri[i];
48  m_kin[i].i = (ii ? ii[i] : 0.0);
49  }
50 
51  if (!inverse) {
52 
53  kiss_fft(m_planf, m_kin, m_kout);
54 
55  for (int i = 0; i < m_n; ++i) {
56  ro[i] = m_kout[i].r;
57  io[i] = m_kout[i].i;
58  }
59 
60  } else {
61 
62  kiss_fft(m_plani, m_kin, m_kout);
63 
64  double scale = 1.0 / m_n;
65 
66  for (int i = 0; i < m_n; ++i) {
67  ro[i] = m_kout[i].r * scale;
68  io[i] = m_kout[i].i * scale;
69  }
70  }
71  }
72 
73 private:
74  int m_n;
75  kiss_fft_cfg m_planf;
76  kiss_fft_cfg m_plani;
77  kiss_fft_cpx *m_kin;
78  kiss_fft_cpx *m_kout;
79 };
80 
81 FFT::FFT(int n) :
82  m_d(new D(n))
83 {
84 }
85 
87 {
88  delete m_d;
89 }
90 
91 void
92 FFT::process(bool inverse,
93  const double *p_lpRealIn, const double *p_lpImagIn,
94  double *p_lpRealOut, double *p_lpImagOut)
95 {
96  m_d->process(inverse,
97  p_lpRealIn, p_lpImagIn,
98  p_lpRealOut, p_lpImagOut);
99 }
100 
102 {
103 public:
104  D(int n) : m_n(n) {
105  if (n % 2) {
106  throw std::invalid_argument
107  ("nsamples must be even in FFTReal constructor");
108  }
109  m_planf = kiss_fftr_alloc(m_n, 0, NULL, NULL);
110  m_plani = kiss_fftr_alloc(m_n, 1, NULL, NULL);
111  m_c = new kiss_fft_cpx[m_n];
112  }
113 
114  ~D() {
115  kiss_fftr_free(m_planf);
116  kiss_fftr_free(m_plani);
117  delete[] m_c;
118  }
119 
120  void forward(const double *ri, double *ro, double *io) {
121 
122  kiss_fftr(m_planf, ri, m_c);
123 
124  for (int i = 0; i <= m_n/2; ++i) {
125  ro[i] = m_c[i].r;
126  io[i] = m_c[i].i;
127  }
128 
129  for (int i = 0; i + 1 < m_n/2; ++i) {
130  ro[m_n - i - 1] = ro[i + 1];
131  io[m_n - i - 1] = -io[i + 1];
132  }
133  }
134 
135  void forwardMagnitude(const double *ri, double *mo) {
136 
137  double *io = new double[m_n];
138 
139  forward(ri, mo, io);
140 
141  for (int i = 0; i < m_n; ++i) {
142  mo[i] = sqrt(mo[i] * mo[i] + io[i] * io[i]);
143  }
144 
145  delete[] io;
146  }
147 
148  void inverse(const double *ri, const double *ii, double *ro) {
149 
150  // kiss_fftr.h says
151  // "input freqdata has nfft/2+1 complex points"
152 
153  for (int i = 0; i < m_n/2 + 1; ++i) {
154  m_c[i].r = ri[i];
155  m_c[i].i = ii[i];
156  }
157 
158  kiss_fftri(m_plani, m_c, ro);
159 
160  double scale = 1.0 / m_n;
161 
162  for (int i = 0; i < m_n; ++i) {
163  ro[i] *= scale;
164  }
165  }
166 
167 private:
168  int m_n;
169  kiss_fftr_cfg m_planf;
170  kiss_fftr_cfg m_plani;
171  kiss_fft_cpx *m_c;
172 };
173 
175  m_d(new D(n))
176 {
177 }
178 
180 {
181  delete m_d;
182 }
183 
184 void
185 FFTReal::forward(const double *ri, double *ro, double *io)
186 {
187  m_d->forward(ri, ro, io);
188 }
189 
190 void
191 FFTReal::forwardMagnitude(const double *ri, double *mo)
192 {
193  m_d->forwardMagnitude(ri, mo);
194 }
195 
196 void
197 FFTReal::inverse(const double *ri, const double *ii, double *ro)
198 {
199  m_d->inverse(ri, ii, ro);
200 }
201 
202 
203 
void forwardMagnitude(const double *ri, double *mo)
Definition: FFT.cpp:135
D * m_d
Definition: FFT.h:101
void inverse(const double *ri, const double *ii, double *ro)
Definition: FFT.cpp:148
~D()
Definition: FFT.cpp:33
~FFTReal()
Definition: FFT.cpp:179
Definition: FFT.cpp:23
kiss_fft_cpx * m_kout
Definition: FFT.cpp:78
void process(bool inverse, const double *ri, const double *ii, double *ro, double *io)
Definition: FFT.cpp:40
FFTReal(int nsamples)
Construct an FFT object to carry out real-to-complex transforms of size nsamples. ...
Definition: FFT.cpp:174
kiss_fft_cpx * m_c
Definition: FFT.cpp:171
void inverse(const double *realIn, const double *imagIn, double *realOut)
Carry out an inverse real transform (i.e.
Definition: FFT.cpp:197
D * m_d
Definition: FFT.h:42
void forward(const double *realIn, double *realOut, double *imagOut)
Carry out a forward real-to-complex transform of size nsamples, where nsamples is the value provided ...
Definition: FFT.cpp:185
#define NULL
Definition: Filter.h:20
void process(bool inverse, const double *realIn, const double *imagIn, double *realOut, double *imagOut)
Carry out a forward or inverse transform (depending on the value of inverse) of size nsamples...
Definition: FFT.cpp:92
kiss_fftr_cfg m_plani
Definition: FFT.cpp:170
kiss_fftr_cfg m_planf
Definition: FFT.cpp:169
void forward(const double *ri, double *ro, double *io)
Definition: FFT.cpp:120
kiss_fft_cfg m_planf
Definition: FFT.cpp:75
kiss_fft_cfg m_plani
Definition: FFT.cpp:76
int m_n
Definition: FFT.cpp:74
D(int n)
Definition: FFT.cpp:26
D(int n)
Definition: FFT.cpp:104
~FFT()
Definition: FFT.cpp:86
FFT(int nsamples)
Construct an FFT object to carry out complex-to-complex transforms of size nsamples.
Definition: FFT.cpp:81
int m_n
Definition: FFT.cpp:168
kiss_fft_cpx * m_kin
Definition: FFT.cpp:77
void forwardMagnitude(const double *realIn, double *magOut)
Carry out a forward real-to-complex transform of size nsamples, where nsamples is the value provided ...
Definition: FFT.cpp:191