FFT problem: FFT -> some manipulation -> IFFT results in one-element shift

46 views (last 30 days)
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 Comments
Oshani
Oshani about 9 hours ago
Moved: Matt J about 1 hour ago

My random noise is always even because we define our random noise in x- minute chunks (sorry that's unclear but what I mean is it's always eg 10 or 15 minutes long, never 10.5 or anything like that, and the update rate is always 1/7th or 1/8th of a second). I will change the code to make it more generalizable though just in case we decide to switch it up in the future!

Oshani
Oshani 17 minutes ago
update: changed my frequency definition and moved it to before any processing takes place which unfortunately did not fix the problem.

Sign in to comment.

Answers (0)

Products


Release

R2025a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!