// // Programmer: Craig Stuart Sapp // Creation Date: Sun Jul 8 14:50:17 PDT 2001 // Last Modified: Sun Jul 8 14:50:20 PDT 2001 // Filename: ...sig/code/utilities/transforms/transforms-templates.cpp // Syntax: C++ // // Description: Basic signal <-> frequency transforms. // #ifndef _TRANSFORMS_TEMPLATES_CPP_INCLUDED #define _TRANSFORMS_TEMPLATES_CPP_INCLUDED #include "Block.h" #include "ComplexD.h" #include #include #include #ifndef OLDCPP #include using namespace std; #else #include #endif #include "transforms.h" #include "transforms-private.h" ////////////////////////////// // // RealFFT -- Calculates the complex spectrum of a real signal. // Uses a N/2 length FFT and therfore more efficient // than the plain FFT for purely real signals. Can accept // Block or Block. Input signal length must // be a power of 2. // template Block& RealFFT(Block& outputSpectrum, const Block& realSignal) { int N = realSignal.getSize()/2; if (!poweroftwo(N)) { cerr << "You can only take the RealFFT of a block with length being" << " a power of 2.\nRequested transform length: " << N << endl; exit(1); } const Block& x = realSignal; outputSpectrum.setSize(2*N); Block& X = outputSpectrum; type pi = 4.0 * atan(1.0); int n, k; // Pre-twiddle: Interpret x as N interleaved complex samples Block y(N); for (n=0; n Y(N); FFT(Y, y); // Post-twiddle: for (k=0; k and Block as input. // template Block& DirectMDCT(Block& outputSpectrum, const Block& inputSignal) { int N = inputSignal.getSize(); Block& x = inputSignal; outputSpectrum.setSize(N/2); Block& X = outputSpectrum; type n0 = ((type)N / 2.0 + 1.0) / 2.0; type tpn = 2.0 * (4.0 * atan(1.0)) / N; type sum; int k, n; if (N%2 != 0) { cerr << "Error: must have an even number of elements to do an MDCT." << endl; exit(1); } for (k=0; k and Block as input. // template Block& FastMDCT(Block& outputSpectrum, const Block& inputSignal) { Block& x = inputSignal; int N = inputSignal.getSize(); outputSpectrum.setSize(N/2); Block& X = outputSpectrum; type pi = 4.0 * atan(1.0); type tpn0 = 2.0 * pi * (N/2.0+1.0)/2.0; int n, k; // Pre-twiddle $$y(n) = x(n)\cos(\pi n/N) - j x(n)\sin(\pi n/N)$$ Block y(N); for (n=0; n Y; FFT(Y, y); // only going to need first half of spectrum // Post-twiddle $$X(k) = \re{Y(k)}\cos(2\pi n_0(k+1/2)/N) + // \im{Y(k)}\sin(2\pi n_0(k+1/2)/N) for (k=0; k Block& MDCT(Block& outputSpectrum, const Block& inputSignal) { int N = inputSignal.getSize(); if (poweroftwo(N) && N > 127) { return FastMDCT(outputSpectrum, inputSignal); } else { return DirectMDCT(outputSpectrum, inputSignal); } } ////////////////////////////// // // DirectIMDCT -- Inverse Modified Discrete Cosine Transform O(N^2). // Can handle Block and Block as input. // template Block& DirectIMDCT(Block& outputSignal, const Block& inputSpectrum) { int N = inputSpectrum.getSize()*2; Block& X = inputSpectrum; outputSignal.setSize(N); Block& x = outputSignal; type n0 = ((type)N/2.0+1.0)/2.0; type pi = 4.0 * atan(1.0); type tpn = 2.0 * pi/N; int k, n; type sum; for (n=0; n and Block as input. // template Block& FastIMDCT(Block& outputSignal, const Block& inputSpectrum) { int N = inputSpectrum.getSize()*2; Block fullSpectrum(N); Block& X = fullSpectrum; outputSignal.setSize(N); Block& x = outputSignal; type n0 = ((type)N/2.0+1.0)/2.0; type pn = (4.0 * atan(1.0))/N; type tpn0 = 2.0 * pn * n0; int n, k; // recreate full spectrum for (n=0; n Z(N); for (k=0; k z; IFFT(z, Z); // Post-twiddle:$$x(n) = 2 \re{z(n)}\cos(2\pi(n+n_0)/2N) - // 2 \im{z(n)}\sin(2\pi(n+n_0)/2N)  for (n=0; n Block& IMDCT(Block& outputSignal, const Block& inputSpectrum) { int N = inputSpectrum.getSize(); if (poweroftwo(N) && N > 127) { return FastIMDCT(outputSignal, inputSpectrum); } else { return DirectIMDCT(outputSignal, inputSpectrum); } } #endif /* _TRANSFORMS_TEMPLATES_CPP_INCLUDED */ // md5sum: bbf1e3eac2b687cd841ea1f871a833a8 transforms-templates.cpp [20050403]