C ++ Audio Processing - FFT

Question

I'm probably going to ask it wrong and make myself look very stupid, but here goes:

I am trying to do some manipulation and sound processing in a wav file. Now I can read all the data (including the header), but I need the data to be in frequency, and for this I need to use the FFT.

I searched the Internet high and low and found it, and the example was taken out of the book “Numerical Recipes in C”, however I changed it to use vectors instead of arrays. So here is the problem:

I was given (as an example to use) a series of numbers and a sampling rate:

X = {50, 206, -100, -65, -50, -6, 100, -135} 

Sampling Rate: 8000 Sample Number: 8

And so I have to answer this:

  0Hz A=0 D=1.57079633 1000Hz A=50 D=1.57079633 2000HZ A=100 D=0 3000HZ A=100 D=0 4000HZ A=0 D=3.14159265 

However, the code that I rewrote compiles, trying to introduce these numbers into the equation (function), I get a segmentation error. Is there something wrong with my code or is the sampling rate too high? (The algorithm does not segment using a much lower sampling rate). Here is the code:

 #include <iostream> #include <math.h> #include <vector> using namespace std; #define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr; #define pi 3.14159 void ComplexFFT(vector<float> &realData, vector<float> &actualData, unsigned long sample_num, unsigned int sample_rate, int sign) { unsigned long n, mmax, m, j, istep, i; double wtemp,wr,wpr,wpi,wi,theta,tempr,tempi; // CHECK TO SEE IF VECTOR IS EMPTY; actualData.resize(2*sample_rate, 0); for(n=0; (n < sample_rate); n++) { if(n < sample_num) { actualData[2*n] = realData[n]; }else{ actualData[2*n] = 0; actualData[2*n+1] = 0; } } // Binary Inversion n = sample_rate << 1; j = 0; for(i=0; (i< n /2); i+=2) { if(j > i) { SWAP(actualData[j], actualData[i]); SWAP(actualData[j+1], actualData[i+1]); if((j/2)<(n/4)) { SWAP(actualData[(n-(i+2))], actualData[(n-(j+2))]); SWAP(actualData[(n-(i+2))+1], actualData[(n-(j+2))+1]); } } m = n >> 1; while (m >= 2 && j >= m) { j -= m; m >>= 1; } j += m; } mmax=2; while(n > mmax) { istep = mmax << 1; theta = sign * (2*pi/mmax); wtemp = sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi = sin(theta); wr = 1.0; wi = 0.0; for(m=1; (m < mmax); m+=2) { for(i=m; (i <= n); i += istep) { j = i*mmax; tempr = wr*actualData[j-1]-wi*actualData[j]; tempi = wr*actualData[j]+wi*actualData[j-1]; actualData[j-1] = actualData[i-1] - tempr; actualData[j] = actualData[i]-tempi; actualData[i-1] += tempr; actualData[i] += tempi; } wr = (wtemp=wr)*wpr-wi*wpi+wr; wi = wi*wpr+wtemp*wpi+wi; } mmax = istep; } // determine if the fundamental frequency int fundemental_frequency = 0; for(i=2; (i <= sample_rate); i+=2) { if((pow(actualData[i], 2)+pow(actualData[i+1], 2)) > pow(actualData[fundemental_frequency], 2)+pow(actualData[fundemental_frequency+1], 2)) { fundemental_frequency = i; } } } int main(int argc, char *argv[]) { vector<float> numbers; vector<float> realNumbers; numbers.push_back(50); numbers.push_back(206); numbers.push_back(-100); numbers.push_back(-65); numbers.push_back(-50); numbers.push_back(-6); numbers.push_back(100); numbers.push_back(-135); ComplexFFT(numbers, realNumbers, 8, 8000, 0); for(int i=0; (i < realNumbers.size()); i++) { cout << realNumbers[i] << "\n"; } } 

Another thing (I know this sounds silly), but I really don’t know what is expected from the "int sign" This is passed through the ComplexFFT function, this is where I can be wrong.

Does anyone have any suggestions or solutions to this problem?

Thanks:)

Answer 1

I think the problem is errors in the way you translated the algorithm.

• Did you mean to initialize j to 1 , not 0 ?

• for(i = 0; (i < n/2); i += 2) should be for (i = 1; i < n; i += 2) .

• Your SWAP should be

   SWAP(actualData[j - 1], actualData[i - 1]); SWAP(actualData[j], actualData[i]); 

• What is SWAP for? I do not think they are needed.

   if((j/2)<(n/4)) { SWAP(actualData[(n-(i+2))], actualData[(n-(j+2))]); SWAP(actualData[(n-(i+2))+1], actualData[(n-(j+2))+1]); } 

• j >= m in while (m >= 2 && j >= m) should be j > m if you intended to make a bit reversal.

• In the code that implements the Danielson-Lanczos section, are you sure that j = i*mmax; should not be an addition, i.e. j = i + mmax; ?

---

In addition, you can do a lot to simplify the code.

Using a SWAP macro should be discouraged when you can just use std::swap ... I was going to suggest std::swap_ranges , but then I realized that you only need to exchange the real parts, since your data is all real (your imaginary parts of the time series - all 0 ):

 std::swap(actualData[j - 1], actualData[i - 1]); 

You can simplify the whole thing using std::complex .

Answer 2

I expect this to be the resize of your vector.

One possibility: maybe resizing will create temporary objects on the stack before moving them back to the heap, I think.

Answer 3

FFT in numerical recipes in C uses the Cooley-Tukey Algorithm , so in response to your question at the end of int sign , the same procedure is passed to calculate both the direct ( sign=-1 ) and the inverse ( sign=1 ) FFT. This is similar to how you use sign when defining theta = sign * (2*pi/mmax) .