Something is wrong with my FFT results

I've been over this for a couple of hours and I can't find where the error is. I just want to confirm that at each step of the process that I'm doing it correctly so all of the outputs should match.

1. Import loudspeaker response data in magnitude and phase.
2. Convert to complex number.
3. Interpolate.
4. IFFT

The first three steps return matching plots, but the plot of the fourth step has something wrong with it. What am I doing wrong??

Link to PRX615M_Main.csv.

And here's the magPhase2complex function I made:

function Z = magPhase2complex(mag_dB,phase_deg)
    Z = [db2mag(mag_dB) .* exp(1j*(deg2rad(phase_deg)))];
end

some n

% Get some Data
CSVmainFile = '/Users/nathanlively/Downloads/PRX615M_Main.csv';
M_TF = readtable(CSVmainFile,'ReadVariableNames',false);
M_TF.Properties.VariableNames = {'frequency','magnitude','phase','coherence'};

% Convert to complex
M_TF.Z = magPhase2complex(M_TF.magnitude,M_TF.phase);

% Interpolate and fill missing
newFs = 192000;
main = timetable('SampleRate',newFs);
main.Hz = linspace(0,newFs,newFs)';
main.Z = interp1(M_TF.frequency,M_TF.Z,main.Hz,'pchip',NaN);
main.Z = fillmissing(main.Z,'nearest','EndValues','extrap');
main.magnitude = mag2db(abs(main.Z));

% IFFT
main.IR = real(ifft(main.Z)) * height(main);
main.Zrecon = fft(main.IR) / height(main); 
main.magnitudeRecon = mag2db(2*abs(main.Zrecon)); 

semilogx(M_TF.frequency,M_TF.magnitude, main.Hz,main.magnitude, main.Hz,main.magnitudeRecon)
xlim([20 20000]);ylim([max(M_TF.magnitude)-10 max(M_TF.magnitude)]);
legend('original','interp1','reconstructed magnitude','FontSize',22)

[Paul]

The code generates an error for the given data set ....

CSVmainFile = 'https://s3.amazonaws.com/appforest_uf/f1621457185394x582849618335596900/JBL%20PRX65M%203.252m.csv';
M_TF = readtable(CSVmainFile,'ReadVariableNames',false);
M_TF.Properties.VariableNames = {'frequency','magnitude','phase','coherence'};
% Convert to complex
M_TF.Z = magPhase2complex(M_TF.magnitude,M_TF.phase);
% Interpolate and fill missing
newFs = 192000;
main = timetable('SampleRate',newFs);
main.Hz = linspace(0,newFs,newFs)';
main.Z = interp1(M_TF.frequency,M_TF.Z,main.Hz,'pchip',NaN);

Error using griddedInterpolant
The sample points must be finite.

Error in interp1 (line 170)
F = griddedInterpolant(X,V(:,1),method);

main.Z = fillmissing(main.Z,'nearest','EndValues','extrap');
main.magnitude = mag2db(abs(main.Z));


function Z = magPhase2complex(mag_dB,phase_deg)
    Z = [db2mag(mag_dB) .* exp(1j*(deg2rad(phase_deg)))];
end

[Nathan Lively]

Thanks so much Paul. I re-ran your code in MATLAB and not only does it not generate the error you are seeing, but now the plots all match. Argh!! ¯\_(ツ)_/¯

Maybe it was a case of not clearing some variable. Unsure.

Response to Paul video

Here's the exact code I used:

% Get some Data
CSVmainFile = 'https://s3.amazonaws.com/appforest_uf/f1621457185394x582849618335596900/JBL%20PRX65M%203.252m.csv';
M_TF = readtable(CSVmainFile,'ReadVariableNames',false);
M_TF.Properties.VariableNames = {'frequency','magnitude','phase','coherence'};

% Convert to complex
M_TF.Z = magPhase2complex(M_TF.magnitude,M_TF.phase);

% Interpolate and fill missing
newFs = 192000;
main = timetable('SampleRate',newFs);
main.Hz = linspace(0,newFs,newFs)';
main.Z = interp1(M_TF.frequency,M_TF.Z,main.Hz,'pchip',NaN);
main.Z = fillmissing(main.Z,'nearest','EndValues','extrap');
main.magnitude = mag2db(abs(main.Z));

% IFFT
main.IR = real(ifft(main.Z)) * height(main);
main.Zrecon = fft(main.IR) / height(main); 
main.magnitudeRecon = mag2db(2*abs(main.Zrecon)); 
semilogx(M_TF.frequency,M_TF.magnitude, main.Hz,main.magnitude, main.Hz,main.magnitudeRecon)
xlim([20 20000]);ylim([max(M_TF.magnitude)-10 max(M_TF.magnitude)]);
legend('original','interp1','reconstructed magnitude','FontSize',22)

[Nathan Lively]

Well, never mind. I opened this up again just now and the problem is back. Should I submit a bug report?

clc;close;clear;

% Get some Data
CSVmainFile = 'https://s3.amazonaws.com/appforest_uf/f1637509399341x777928112764030700/PS15R2%20WB.csv';
M_TF = readtable(CSVmainFile,'ReadVariableNames',false);
M_TF.Properties.VariableNames = {'frequency','magnitude','phase','coherence'};

% Convert to complex
M_TF.Z = magPhase2complex(M_TF.magnitude,M_TF.phase);

% Interpolate and fill missing
newFs = 192000;
main = timetable('SampleRate',newFs);
main.Hz = linspace(0,newFs,newFs)';
main.Z = interp1(M_TF.frequency,M_TF.Z,main.Hz,'pchip',NaN);
main.Z = fillmissing(main.Z,'nearest','EndValues','extrap');
main.coherence = interp1(M_TF.frequency,M_TF.coherence,main.Hz);
main.coherence(isnan(main.coherence)) = 0;
main.magnitude = mag2db(abs(main.Z));

% IFFT
main.IR = real(ifft(main.Z)) * height(main);
main.Zrecon = fft(main.IR) / height(main); 
main.magnitudeRecon = mag2db(2*abs(main.Zrecon)); 

semilogx(M_TF.frequency,M_TF.magnitude, main.Hz,main.magnitude, main.Hz,main.magnitudeRecon)
xlim([20 20000]);ylim([max(M_TF.magnitude)-30 max(M_TF.magnitude)]);
legend('original','interp1','reconstructed magnitude','FontSize',22)

[Nathan Lively]

Maybe the FFT function doesn't like the interpolated values?

[Nathan Lively]

@Paul I received an email with another comment from you, but I don't see it here. Anway, I'll reply to it.

Does the orginal data set contain more or less than 192000 points?

No, this is one of the problems with the original data set. Sometimes it does not start at zero, so I'm trying to fix that, upsample, and standardize the frequency spacing (1Hz) so that everything coming in get cleaned up and matched.

And is the original data set equally spaced in frequency?

No. This is another reason why I'm using interpolation.

I suspect that this line should be changed to get a unit freuqency spacing, though probably not a noticeable difference .

%main.Hz = linspace(0,newFs,newFs )';
main.Hz = linspace(0,newFs-1,newFs )';

Yes! You are totally right. But that does not fix the problem, unfortunately.

[Chris Turnes]

Perhaps you want to do

main.IR = ifft(main.Z, 'symmetric') * height(main);

and

main.magnitudeRecon = mag2db(abs(main.Zrecon)); % No factor of 2

This appears to give matching amplitudes (at least, for me). The 'symmetric' part is because you're not interpolating the portion of the spectrum that corresponds to negative frequencies (you're just putting in zeros there) and the removal of the 2 is because I can't follow why it was there to begin with

[Nathan Lively]

Maybe I found a solution?

Seems that when I change:

main.Z = interp1(M_TF.frequency,M_TF.Z,main.Hz,'pchip',NaN);

to

main.Z = interp1(M_TF.frequency,M_TF.Z,main.Hz,'pchip',0);

Then the FFT function works better. I wish I knew if this was the best way to do it, but interp1 is the only way I have found that works when you have missing on both ends.

[Nathan Lively]

Thanks @Chris Turnes! I'm not sure why the 2 goes there either. I'm taking a class on FFT and it says, "Multiply by 2 because amplitude is split between positive and negative frequencies."

Your suggestion works, though.

[Paul]

@Chris Turnes

The comment about using the 'symmetric' option suggests that if that option is used then the "upper half" of the input to ifft() doesn't matter. Am I interpreting that comment correctly? If not, can you clarify? If so, I don't think that's the case and it's not stated as such in the doc.

@Nathan Lively

I think I'm still not clear on what the goal is. As I understand it, you have data of a spectrum in the frequency domain that is not equally spaced, doesn't start at f = 0, and likely doesn't exactly end at Fs/2, which looks like maybe it should be 48000 Hz. If that's the data you have, what is the goal? Generate a new spectrum, that starts at f =,0, is equally spaced at 1 Hz, and finishes at Fs/2 = 192000/2? If so what is the use of that end product? I'm a bit confused because of the ifft -> fft operation in the code in the question.

Also, the first two rows in the table imported from the .csv file are NaN for magnitude and phase, which is why I was getting the error in the comment above.

[Chris Turnes]

@Paul - I think the documentation states this though, perhaps not as explicitly as you're hinting. You are correct though -- if you specify the 'symmetric' flag, the input is treated as if it is conjugate symmetric so effectively what is in the back half is ignored. This is really meant for cases where, through various operations, you started off with something that should be the transform of a real vector and somehow became close to but not exactly conjugate symmetric. The documentation page has a sentence that says "For example, ifft(Y,'symmetric') treats Y as conjugate symmetric" which doesn't explicitly say that the second half is ignored but implies it.

You can convince yourself this is what is indeed happening fairly easily:

x = randn(1,11);
y = fft(x);
% Fill the second half with stuff that would obviously distort the transform
y(7:end) = realmax;
ifft(y,'symmetric')

ans = 1×11

-0.2682 0.5258 -0.5491 0.7655 -1.5322 -0.5287 0.0856 0.4099 -0.0580 -0.4629 0.7874

norm(ifft(y,'symmetric') - x)

ans = 3.4360e-16

If any of those realmax's were actually used, you'd see non-finites in the results.

[Paul]

Sure enough. I attempted this experiment (different data of course) and didn't get this result. I must have done something wrong. Thanks for pointing this out.

As for the doc, I think it should just state explicitly that the upper half is ignored.

[Chris Turnes]

That's a reasonable request; I'll pass it on to our documentation team.

[Bjorn Gustavsson]

@Chris Turnes: Is it obvious that one should use the first half of a "should-be symmetric" FT that has "gone bad"? Couldn't the numerical errors be rather evenly distributed between the positive and "negative" frequency-components?

[Chris Turnes]

@Bjorn Gustavsson - Sure, that's entirely possible. I'm sure there are arguments for using a different convention than just taking the first half; the option is really meant as a convenience for situations where there is very small round-off occurring, in which case it probably doesn't practically matter what convention is used.

The motivation for the choice of how this is handled very likely came from the fact that MATLAB uses FFTW to compute FFTs, and FFTW's "C2R" algorithm for transforming complex data back into real data follows this exact convention.

[Bjorn Gustavsson]

@Chris Turnes - fair points, the question just struck me as I read your description, and thought that "the optimal procedure surely ought to be the averages of the real and the conjugated complex parts". But I see the use-case and that my flicker of an idea doesn't really apply...

[Nathan Lively]

@Paul Generate a new spectrum, that starts at f =,0, is equally spaced at 1 Hz, and finishes at Fs/2 = 192000/2?

Exactly!

If so what is the use of that end product?

One of the problems I was running into is that I want to window the IFFT data at a specific lenght in milliseconds. With the "dirty" data I am not able to calculate the time vector to go along with the IFFT.

Answer 1

[Nathan Lively]

Although I'm not 100% this is correct, so far the answer I mentioned in this comment seems to be working. I just needed to replace the beginning and end values with zeros.