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