Something is wrong with my FFT results

26 次查看（过去 30 天）

Nathan Lively 2022-3-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1667124-something-is-wrong-with-my-fft-results

回答： Nathan Lively 2022-3-10

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.

Import loudspeaker response data in magnitude and phase.
Convert to complex number.
Interpolate.
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)

16 个评论
显示 14更早的评论隐藏 14更早的评论

Nathan Lively 2022-3-9

编辑：Nathan Lively 2022-3-9

在 MATLAB Online 中打开

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 2022-3-9

在 MATLAB Online 中打开

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)

Bjorn Gustavsson 2022-3-10

@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 2022-3-10

@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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Nathan Lively 2022-3-10

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1667124-something-is-wrong-with-my-fft-results#answer_914779

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.