← Back to branch summary

~ubuntu-branches/ubuntu/utopic/ardour3/utopic

« back to all changes in this revision

Viewing changes to libs/qm-dsp/dsp/transforms/FFT.cpp

  • Committer: Package Import Robot
  • Author(s): Felipe Sateler
  • Date: 2013-09-21 19:05:02 UTC
  • Revision ID: package-import@ubuntu.com-20130921190502-8gsftrku6jnzhd7v
Tags: upstream-3.4~dfsg
Import upstream version 3.4~dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
libs/qm-dsp/dsp/transforms/FFT.cpp
 
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 <cmath>
 
15
 
 
16
#include <iostream>
 
17
 
 
18
FFT::FFT(unsigned int n) :
 
19
    m_n(n),
 
20
    m_private(0)
 
21
{
 
22
    if( !MathUtilities::isPowerOfTwo(m_n) )
 
23
    {
 
24
        std::cerr << "ERROR: FFT: Non-power-of-two FFT size "
 
25
                  << m_n << " not supported in this implementation"
 
26
                  << std::endl;
 
27
        return;
 
28
    }
 
29
}
 
30
 
 
31
FFT::~FFT()
 
32
{
 
33
 
 
34
}
 
35
 
 
36
FFTReal::FFTReal(unsigned int n) :
 
37
    m_n(n),
 
38
    m_private_real(0)
 
39
{
 
40
    m_private_real = new FFT(m_n);
 
41
}
 
42
 
 
43
FFTReal::~FFTReal()
 
44
{
 
45
    delete (FFT *)m_private_real;
 
46
}
 
47
 
 
48
void
 
49
FFTReal::process(bool inverse,
 
50
                 const double *realIn,
 
51
                 double *realOut, double *imagOut)
 
52
{
 
53
    ((FFT *)m_private_real)->process(inverse, realIn, 0, realOut, imagOut);
 
54
}
 
55
 
 
56
static unsigned int numberOfBitsNeeded(unsigned int p_nSamples)
 
57
{       
 
58
    int i;
 
59
 
 
60
    if( p_nSamples < 2 )
 
61
    {
 
62
        return 0;
 
63
    }
 
64
 
 
65
    for ( i=0; ; i++ )
 
66
    {
 
67
        if( p_nSamples & (1 << i) ) return i;
 
68
    }
 
69
}
 
70
 
 
71
static unsigned int reverseBits(unsigned int p_nIndex, unsigned int p_nBits)
 
72
{
 
73
    unsigned int i, rev;
 
74
 
 
75
    for(i=rev=0; i < p_nBits; i++)
 
76
    {
 
77
        rev = (rev << 1) | (p_nIndex & 1);
 
78
        p_nIndex >>= 1;
 
79
    }
 
80
 
 
81
    return rev;
 
82
}
 
83
 
 
84
void
 
85
FFT::process(bool p_bInverseTransform,
 
86
             const double *p_lpRealIn, const double *p_lpImagIn,
 
87
             double *p_lpRealOut, double *p_lpImagOut)
 
88
{
 
89
    if (!p_lpRealIn || !p_lpRealOut || !p_lpImagOut) return;
 
90
 
 
91
//    std::cerr << "FFT::process(" << m_n << "," << p_bInverseTransform << ")" << std::endl;
 
92
 
 
93
    unsigned int NumBits;
 
94
    unsigned int i, j, k, n;
 
95
    unsigned int BlockSize, BlockEnd;
 
96
 
 
97
    double angle_numerator = 2.0 * M_PI;
 
98
    double tr, ti;
 
99
 
 
100
    if( !MathUtilities::isPowerOfTwo(m_n) )
 
101
    {
 
102
        std::cerr << "ERROR: FFT::process: Non-power-of-two FFT size "
 
103
                  << m_n << " not supported in this implementation"
 
104
                  << std::endl;
 
105
        return;
 
106
    }
 
107
 
 
108
    if( p_bInverseTransform ) angle_numerator = -angle_numerator;
 
109
 
 
110
    NumBits = numberOfBitsNeeded ( m_n );
 
111
 
 
112
 
 
113
    for( i=0; i < m_n; i++ )
 
114
    {
 
115
        j = reverseBits ( i, NumBits );
 
116
        p_lpRealOut[j] = p_lpRealIn[i];
 
117
        p_lpImagOut[j] = (p_lpImagIn == 0) ? 0.0 : p_lpImagIn[i];
 
118
    }
 
119
 
 
120
 
 
121
    BlockEnd = 1;
 
122
    for( BlockSize = 2; BlockSize <= m_n; BlockSize <<= 1 )
 
123
    {
 
124
        double delta_angle = angle_numerator / (double)BlockSize;
 
125
        double sm2 = -sin ( -2 * delta_angle );
 
126
        double sm1 = -sin ( -delta_angle );
 
127
        double cm2 = cos ( -2 * delta_angle );
 
128
        double cm1 = cos ( -delta_angle );
 
129
        double w = 2 * cm1;
 
130
        double ar[3], ai[3];
 
131
 
 
132
        for( i=0; i < m_n; i += BlockSize )
 
133
        {
 
134
 
 
135
            ar[2] = cm2;
 
136
            ar[1] = cm1;
 
137
 
 
138
            ai[2] = sm2;
 
139
            ai[1] = sm1;
 
140
 
 
141
            for ( j=i, n=0; n < BlockEnd; j++, n++ )
 
142
            {
 
143
 
 
144
                ar[0] = w*ar[1] - ar[2];
 
145
                ar[2] = ar[1];
 
146
                ar[1] = ar[0];
 
147
 
 
148
                ai[0] = w*ai[1] - ai[2];
 
149
                ai[2] = ai[1];
 
150
                ai[1] = ai[0];
 
151
 
 
152
                k = j + BlockEnd;
 
153
                tr = ar[0]*p_lpRealOut[k] - ai[0]*p_lpImagOut[k];
 
154
                ti = ar[0]*p_lpImagOut[k] + ai[0]*p_lpRealOut[k];
 
155
 
 
156
                p_lpRealOut[k] = p_lpRealOut[j] - tr;
 
157
                p_lpImagOut[k] = p_lpImagOut[j] - ti;
 
158
 
 
159
                p_lpRealOut[j] += tr;
 
160
                p_lpImagOut[j] += ti;
 
161
 
 
162
            }
 
163
        }
 
164
 
 
165
        BlockEnd = BlockSize;
 
166
 
 
167
    }
 
168
 
 
169
 
 
170
    if( p_bInverseTransform )
 
171
    {
 
172
        double denom = (double)m_n;
 
173
 
 
174
        for ( i=0; i < m_n; i++ )
 
175
        {
 
176
            p_lpRealOut[i] /= denom;
 
177
            p_lpImagOut[i] /= denom;
 
178
        }
 
179
    }
 
180
}
 
181