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!!