FFTW result different for one vs two dimensional matrix

41 Views Asked by At

I am trying to calculate the Fourier transform of each row of a matrix, and I am getting different outputs depending on whether the input matrix has one or two rows. I am not able to just take the 2D Fourier transform directly, as the result of the modification of the first row will affect the second once I put in other code that I need, so I need to take the Fourier transform sequentially. Right now, I am just trying to return the input after taking the Fourier transform twice. Here is my code:

std::vector<std::vector<std::complex<double>>> output_mat(numRows, std::vector<std::complex<double>>(numCols));

mxArray *outputMatrix = mxCreateDoubleMatrix(numRows, numCols, mxCOMPLEX);
double *realPart = mxGetPr(outputMatrix);
double *imagPart = mxGetPi(outputMatrix);

std::vector<std::complex<double>> u1(numElements_em);
std::vector<std::complex<double>> u1_pre(numElements_em);
std::vector<std::complex<double>> uip(numCols);
std::vector<std::complex<double>> uhalf_loop(numCols);
std::vector<std::complex<double>> uhalf(numCols);

// ~~~~~~~~~~~~~~~~~~~ BEGINNING LOOP! ~~~~~~~~~~~~~~~~~~~~~~

for (int pulse = 0; pulse < num_pulses; ++pulse) {

    fftw_complex* in_fft = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * numCols);
    fftw_complex* out_fft = (fftw_complex*)fftw_malloc(sizeof(fftw_complex) * numCols);

    fftw_plan forward_plan = fftw_plan_dft_1d(numCols, in_fft, out_fft, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_plan backward_plan = fftw_plan_dft_1d(numCols, out_fft, in_fft, FFTW_BACKWARD, FFTW_ESTIMATE);

    u1 = usArray_mat[pulse];

    std::cout << pulse << std::endl;
    std::cout << u1.size() << std::endl;

    //u1_pre = u1;

    // Initialize in_fft with the real data and set imaginary part to 0
    for (int i = 0; i < u1.size(); i++) {
        in_fft[i][0] = std::real(u1[i]);
        in_fft[i][1] = std::imag(u1[i]); // Set imaginary part to 0
    }

    fftw_execute(forward_plan);

    for (int i = 0; i < u1.size(); ++i) {
        uip[i] = std::complex<double>(out_fft[i][0], out_fft[i][1]);
        uhalf_loop[i] = uip[i]; 
    }

    for (int i = 0; i < uhalf_loop.size(); i++) {
        out_fft[i][0] = std::real(uhalf_loop[i])/uhalf_loop.size();
        out_fft[i][1] = std::imag(uhalf_loop[i])/uhalf_loop.size(); // Set imaginary part to 0
    }

    fftw_execute(backward_plan);

    for (int i = 0; i < u1.size(); ++i) {
        uhalf[i] = std::complex<double>(in_fft[i][0], in_fft[i][1]);
    }

    output_mat[pulse] = uhalf;

    fftw_destroy_plan(forward_plan);
    fftw_destroy_plan(backward_plan);
    fftw_free(in_fft);
    fftw_free(out_fft);
 
}


// ~~~~~~~~~~~~~~~~~~~ Creating output matrices ~~~~~~~~~~~~~~~~~~~~~~

// Create output mxArray

for (mwIndex i = 0; i < numRows; i++) { 
    for (mwIndex j = 0; j < numCols; j++) {
        //std::complex<double> complexValue = us_vector;
        realPart[i + numRows * j] = std::real(output_mat[i][j]);
        imagPart[i + numRows * j] = std::imag(output_mat[i][j]);;
    }
}

When I tried the above for a 1D case (i.e. one row x 32768 cols representing one pulse), I got the input back as expected. But when I tried the 2D case (i.e. two rows x 32768 cols representing two pulses), the spectrum had changed and the time domain of the pulse looks like it had been clipped. So, somehow, I am changing the output from the first pulse when I include the second pulse. I have a feeling that I am missing something very simple, but I have tried looking everywhere and have been staring at this code for a while, so I would appreciate any help anyone could give. Thanks!!

0

There are 0 best solutions below