FFT problem: FFT -> some manipulation -> IFFT results in one-element shift
46 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
I am trying to take some random noise (produced by a different script), Fourier transform it, multiply it by a filter, and then do an inverse Fourier transform to get some random noise data back. When I was testing my script with a filter which has a value of 1 for all frequencies (i.e. should give back the original data after IFFT), I realized all of the resulting values are shifted by 1 element. I.e., if I plot the original vs multiplied ('calibrated') noise it looks like this:
but then if I just circshift one of the arrays, either the noise given or the post-FT noise, by 1 (in opposite directions), it looks like this:

(and the values of the two arrays are equal to within 1e-15 tolerance, which is good enough for my purpose). I suspect the problem is how I'm setting up my frequency arrays/interpolating my actual function I need to multiply by, but I'm pretty much stuck after that--any suggestions would be appreciated.
More in-depth info
Full script is attached (spaghetti code beware, apologies in advance), but it has a lot of plots/troubleshooting, so here is the relevant code:
%%read files
[tempFile, tempPath] = uigetfile('*.csv', 'Select the calibration data file');
calData = readmatrix(fullfile(tempPath, tempFile));
[tempFile2, tempPath2] = uigetfile('*.txt', 'Select the random noise');
randNoise = readmatrix(fullfile(tempPath2, tempFile2));
%% Sort and format calibration daya
%sort calibration data by period
clDatSorted=sortrows(calData,[1]);
%extract data
periods = clDatSorted(:,1);
ampsRequired = clDatSorted(:,2);
%mirror the function in the negative frequency domain due to the nature of
%the FT function
ampsRequired = [flip(ampsRequired);ampsRequired];
periods = [flip(periods)*-1;periods];
freqs = 1./periods;
%% Fourier transform the random noise and plot
fourierNoise = fft(randNoise);
fourierNoise = fftshift(fourierNoise); %fftshift used here to center zero-frequency value
Ts = inputdlg('What is the sample period?'); %sample period; user defined because of the way the script to generate my noise is set up
Ts=str2num(Ts{1}); %convert to a number
f= ((-length(fourierNoise)/2):((length(fourierNoise)/2)-1))*fs/length(fourierNoise);
%% interpolate calibration function
%this is going to be important because my calibration function is manually
%obtained so I will need to extrapolate *and* interpolate, unfortunately--I
%can't just give it a smooth function due to the nature of the data
interpolatedCal=interp1(freqs(2:end-1),ampsRequired(2:end-1),f,'linear','extrap');
%% calibrate
calibratedFourierNoise = interpolatedCal .* fourierNoise';
calibratedNoise = ifft(ifftshift(calibratedFourierNoise));
if size(calibratedNoise)~=size(randNoise)
calibratedNoise = flip(calibratedNoise); %to fix row-column swaps--I know this needs a better fix
end
And
ismembertol(circshift(randNoise,[0, 1]),calibratedNoise,1e-15)
produces 1's everywhere.
I think the line where I define my frequency:
f= ((-length(fourierNoise)/2):((length(fourierNoise)/2)-1))*fs/length(fourierNoise);
is the root of the problem but I've tried playing around with it and keep encountering the same issue.
In particular, if I do
[p,fp] = pspectrum(randNoise);
the maximum frequency it gives me is 3.14 for this data, and I'm assuming that the max frequency for the Fourier transform is 1/2 the number of random noise values per second--I feel like this is an issue with my lack of deep understanding of Fourier theory tbh, but I can't pinpoint the exact issue.
Full script/test calibration file with ones everywhere (i.e. flat filter)/sample noise/sample calibration (with ones everywhere as used for the test in sampleCalOneEverywhere and as the filter I'll aactually be using in sampleCal1) are attached.
(Not super relevant but the final filter will be an empirically-obtained low-pass filter with some added stuff.)
3 Kommentare
Antworten (0)
Siehe auch
Kategorien
Mehr zu Bartlett finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!