FFT result dependent on number of interpolation points
Show older comments
I use interpft to interpolate a time series data (with 4hrs intervals recorded for 96 hours) and use fft to find the frequency/period. However, my results are dependent on the parameter N (number of interpolation points) and the interpolated data is very close but do not pass from the original data points. I am new in fft analysis, appreciate any help/advice.
Accepted Answer
The interpft function is intended to interpolate the Fourier transform, and it does a reasonably good approximation. I checked it against the results of the interp1 function. I believe that what you are seeing is simply an artefact of the plot function, since the innterpolated data appear to exactly match the original data at the same points. There are small differences, however they are in the range of floating-point approximation error (approsimately the same as would be returned by the eps function).
Your code works —
% clc
% clear all
% close all
T = 4 * 3600; % Sampling period [s]- 4hr time steps
Fs = 1/T; % Sampling frequency [1/s] - 4hr time steps
L = 3600*96; % Length of sampling [s]
tdata = 0:T:24*T;
N = 50; % number of interpolation points
dt = (tdata(2)-tdata(1))*numel(tdata)/N;
% dt = (tdata(2)-tdata(1))/N;
tmin = 0;
tmax = L;
tinterp = tmin:dt:tmax;
% tinterp = (0:1000:L-1); % Time vector
% use interpft to perform interpolation in frequency domain after
% interpolation it inverse transforms to frequency domain back to original
% domain
X = [1 0.916546112000000 0.897504521000000 0.955244123000000 0.870361664000000 0.901520531000000 0.970887620000000 0.915627954000000 0.884539850000000 0.976460857000000 0.872153443000000 0.923043135000000 1.02580986200000 0.941502448000000 0.904196842000000 0.984450007000000 0.859495215000000 0.916077493000000 0.997053985000000 0.912399374000000 0.878544039000000 0.969741146000000 0.885025052000000 0.926594672000000 0.968193226000000];
Xinterp = interpft(X,N);
format longG
Check = [X; Xinterp(1:2:end)];
disp(Check)
CheckDifference = diff(Check);
disp(CheckDifference)
format short
Xinterp = Xinterp(1:numel(tinterp));
Xinterp = Xinterp -mean(Xinterp);
Y = fft(Xinterp);
plot(tdata/3600,X-mean(X),'o')
hold on
plot(tinterp/3600,Xinterp,'.-');
xlabel('Time [hrs]')
ylabel('Normalized TEER')
xlim([0 96]); xticks(0:12:96);
title('Interpolated Time Series Data');
legend('Data','FT Interpolation')
figure
Y_oneside = Y(1:N/2);
f = Fs * (0:N/2-1)/N;
Y_meg = abs(Y_oneside)/(N/2);
plot(f,Y_meg)
xlabel('Frequency [Hz]');
ylabel('Amplitude');
title('Frequency-domain plot')
[M,I] = max(Y_meg) % locate the peak of frequency domain plot
P = 1/f(I); % Period [s]
P_hrs = P/3600 % Period [hrs]
% vq2 = interp1(EXP_1(7,:),v,xq,'spline');
% plot(x,v,'o',xq,vq2,':.');
% % xlim([0 5E-11])
Xsize = numel(X)
% tinterp = linspace(min(tdata), max(tdata), N);
Xinterp = interp1(tdata, X, tinterp)
format longG
Check = [X; Xinterp(1:2:end)];
disp(Check)
CheckDifference = diff(Check);
disp(CheckDifference)
format short
figure
plot(tdata/3600, X-mean(X),'o')
hold on
plot(tinterp/3600, Xinterp-mean(Xinterp), '.--')
hold off
grid
xlabel('t')
ylabel('X')
The interp1 funciton gives a slightly more accurate result than interpft, however the difference is negligible.
.
6 Comments
Ali
on 1 Nov 2024
My pleasure!
You are calculating the fft based on the interpolated values. It took a few minutes to understand what the problem was, however thee essence of it was that you are using the original ‘Fs’ for the interpolated values, when it should be twice that value because you are using ‘tinterp’, so instead, use:
Fs2 = 1/(tinterp(2)-tinterp(1))
and substitute ‘Fs2’ for ‘Fs’ in the ‘f’ calculation, and it works (and agrees with my separate analysis).
The number of interpolation points should not make nuch — if any — significant difference. When I use the function I wrote to perform discrete Fourier transforms (included here), and use it with the original and interpolated data, the Fourier transforms overlap tot a large extent.
% clc
% clear all
% close all
T = 4 * 3600; % Sampling period [s]- 4hr time steps
Fs = 1/T; % Sampling frequency [1/s] - 4hr time steps
L = 3600*96; % Length of sampling [s]
tdata = 0:T:24*T;
N = 50; % number of interpolation points
dt = (tdata(2)-tdata(1))*numel(tdata)/N;
% dt = (tdata(2)-tdata(1))/N;
tmin = 0;
tmax = L;
tinterp = tmin:dt:tmax;
% tinterp = (0:1000:L-1); % Time vector
% use interpft to perform interpolation in frequency domain after
% interpolation it inverse transforms to frequency domain back to original
% domain
X = [1 0.916546112000000 0.897504521000000 0.955244123000000 0.870361664000000 0.901520531000000 0.970887620000000 0.915627954000000 0.884539850000000 0.976460857000000 0.872153443000000 0.923043135000000 1.02580986200000 0.941502448000000 0.904196842000000 0.984450007000000 0.859495215000000 0.916077493000000 0.997053985000000 0.912399374000000 0.878544039000000 0.969741146000000 0.885025052000000 0.926594672000000 0.968193226000000];
Xinterp = interpft(X,N);
format longG
Check = [X; Xinterp(1:2:end)];
disp(Check)
CheckDifference = diff(Check);
disp(CheckDifference)
format short
Xinterp = Xinterp(1:numel(tinterp));
Xinterp = Xinterp -mean(Xinterp);
Y = fft(Xinterp)
plot(tdata/3600,X-mean(X),'o')
hold on
plot(tinterp/3600,Xinterp,'.-');
xlabel('Time [hrs]')
ylabel('Normalized TEER')
xlim([0 96]); xticks(0:12:96);
title('Interpolated Time Series Data');
legend('Data','FT Interpolation')
% Ni = numel(Xinterp)
% Ni = 2*floor(Ni/2)
Fs2 = 1/(tinterp(2)-tinterp(1))
figure
Y_oneside = Y(1:N/2);
f = Fs2 * (0:N/2-1)/N;
Y_meg = abs(Y_oneside)/(N/2);
plot(f,Y_meg)
xlabel('Frequency [Hz]');
ylabel('Amplitude');
title('Frequency-domain plot')
[M,I] = max(Y_meg) % locate the peak of frequency domain plot
P = 1/f(I); % Period [s]
P_hrs = P/3600 % Period [hrs]
% vq2 = interp1(EXP_1(7,:),v,xq,'spline');
% plot(x,v,'o',xq,vq2,':.');
% % xlim([0 5E-11])
% Xsize = numel(X)
% % tinterp = linspace(min(tdata), max(tdata), N);
%
% Xinterp = interp1(tdata, X, tinterp)
%
% format longG
% Check = [X; Xinterp(1:2:end)];
% disp(Check)
% CheckDifference = diff(Check);
% disp(CheckDifference)
% format short
%
% figure
% plot(tdata/3600, X-mean(X),'o')
% hold on
% plot(tinterp/3600, Xinterp-mean(Xinterp), '.--')
% hold off
% grid
% xlabel('t')
% ylabel('X')
[FTs11,Fv1] = FFT1(X,tdata);
[FTs12,Fv2] = FFT1(Xinterp,tinterp);
figure
plot(Fv1, abs(FTs11)*2, '-b', 'DisplayName','Original')
hold on
plot(Fv2, abs(FTs12)*2, '--r', 'DisplayName','Interpolated')
hold off
grid
xlabel('Frequency')
ylabel('Magnitude')
legend('Location','best')
function [FTs1,Fv] = FFT1(s,t)
% Arguments:
% s: Signal Vector Or Matrix
% t: Associated Time Vector
t = t(:);
L = numel(t);
if size(s,2) == L
s = s.';
end
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
Fv = Fs*(0:(NFFT/2))/NFFT;
% Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
Fv = Fv(:);
FTs1 = FTs(Iv,:);
end
.
Ali
on 2 Nov 2024
Thank you Star Strider! Your answer helped me a lot.
Star Strider
on 2 Nov 2024
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.
Ali
on 3 Dec 2024
Star Strider
on 3 Dec 2024
Hello Ali!
To limit the frequencies displayed, you can either use the xlim function, or longical indexing, such as:
fidx = (f >= lower_frequency & (f <= upper_frequency);
figure
plot(f(fidx), Y_meg(fidx))
grid
The parentheses in ‘fidx’ are not striictly necessary, however they make it easier to read.
I do not have acess to the general reference, so I am not certain. You can improve the frequency resolution of the fft by using zero-padding. The best way to do that is to use the nextpow2 function and its result as the ‘nfft’ argument to the fft function. See the documentation for details on how best to use it. The frequencies and the frequency range will not changee, however the frequency resolution will improve. I include a Fourieer transform function that I wrote for my convenience that demonstrates this:
function [FTs1,Fv] = FFT1(s,t)
% Arguments:
% s: Signal Vector Or Matrix
% t: Associated Time Vector
t = t(:);
L = numel(t);
if size(s,2) == L
s = s.';
end
Fs = 1/mean(diff(t));
Fn = Fs/2;
NFFT = 2^nextpow2(L);
FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));
Fv = Fs*(0:(NFFT/2))/NFFT;
% Fv = linspace(0, 1, NFFT/2+1)*Fn;
Iv = 1:numel(Fv);
Fv = Fv(:);
FTs1 = FTs(Iv,:);
end
.
More Answers (2)
Sahas
on 1 Nov 2024
0 votes
As per my understanding, the "interpft" function is used for fourier interpolation by transforming the data into the frequency domain, performing interpolation, and then transforming it back to the time domain. This method does not guarantee that the interpolated curve will pass through the original data points.
If you would like to ensure that the interpolated curve passes through the original points, use MATLAB's "interp1" function for interpolation.
It is also observed that FFT algorithms perform optimally when "N" is a power of two. I would recommend choosing "N" such that it's greater than the original number of data points, which is the length of "X" array in the code provided above.
Refer to the following MathWorks documentation links for more information on the "interpft" and "interp1" functions:
- https://www.mathworks.com/help/matlab/ref/interpft.html
- https://www.mathworks.com/help/matlab/ref/double.interp1.html
I hope this is beneficial!
1 Comment
Ali
on 2 Nov 2024
Thank you Sahas!
Hi Ali,
"interpft(x,n) interpolates the Fourier transform of the function values in X to produce n equally spaced points."
That statement is not really correct insofar as interpft does not interpolate in the frequency domain, certainly not as one would conventionally interpret that statement.
Instead, interpft() performs Frequency Domain Zero Padding (aka FDZP) to interpolate the time domain signal, which is analgous to zero padding in the time domain to interpolate the transform in the frequency domain. A better statement would be
interpft(x,n) zero pads the Discrete Fourier Transform of the function in the the frequency domain and then takes the inverse to produce n equally spaced points in the time domain.
We can illustrate the method as follows.
Define a continuous-time sine wave with a 10 Hz frequency
x = @(t) sin(2*pi*10*t);
Uniform sampling of the signal over 5 seconds with sampling frequencies of 50 and 100 Hz
Fs1 = 50; Ts1 = 1/Fs1;
N1 = 250;
n1 = 0:N1-1;
t1 = n1*Ts1;
x1 = x(t1);
Fs2 = 100; Ts2 = 1/Fs2;
N2 = 500;
n2 = 0:N2-1;
t2 = n2*Ts2;
x2 = x(t2);
Take the DFT of both
X1 = fft(x1);
X2 = fft(x2);
Plot the (amplitude of the) DFTs versus their independent variables.
figure
stem(n1,abs(X1),'DisplayName','X1');
hold on
stem(n2,abs(X2),'DisplayName','X2');
legend
Now, we want to manipulate the data in X1 so that the result matches X2. The way to do that is to add N2 - N1 = 250 zeros between the blue stems, symetrically around the Nyquist point, i.e., zero padding in the frequency domain. The way to think about this is
First get the DFT of X1 centered around 0
X1centered = fftshift(X1);
Now pad with (N2-N1)/2 zeros on both sides (we'd have to be careful here if N2 - N1 were odd ...)
X1padded = [zeros(1,(N2-N1)/2), X1centered, zeros(1,(N2-N1)/2)];
Then shift back
X1padded = ifftshift(X1padded);
Compare X1padded to X2
figure
stem(n2,abs(X1padded),'DisplayName','X1padded')
hold on
stem(n2,abs(X2),'DisplayName','X2')
Now X1padded is aligned with X2, and its amplitude is diffferent by a factor of N1/N2. So we can recover x2 by taking the ifft of X1padded and multiplying the result by N2/N1.
In this case, because N2/N1 is an integer, the resconstructed points from x2recon = N2/N1*ifft(x1padded) will, theoretically, include the exact original points in x1. In practice, the reconstruction won't be perfect because of numerical inaccuracy of ifft(fft()) and roundoff errors in multiplication by N2/N1. If N2/N1 is not an integer, then the original points won't be a subset of the reconstructed points.
In short, interpft* is interpolating the original time domain signal using a sum of complex sinusoids with sampling frequency increased by N2/N1, which would be different than other interpolation techniques. Whether or not the output of interpft provides any additional useful information would depend on how amenable the input data is to iterpolation in the first place and how one uses that output.
*What I've described is essentially how interpft works when the number of requested interpolated points is greater than the number of input points (interpft also takes care to handle the Nyquist point of the DFT correctly when numel(x) is even, which I could ignore for this example because that point is 0). interpft also handles the case where the requested number of interpolated points is fewer than the number input points, which requires some extra processing.
Categories
Find more on Interpolation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
