FFT result dependent on number of interpolation points

32 次查看（过去 30 天）
Ali
Ali 2024-11-1
评论： Star Strider 2024-12-3，16:42
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.

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

采纳的回答

Star Strider
Star Strider 2024-11-1
Ran in:
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)
Columns 1 through 6 1 0.916546112 0.897504521 0.955244123 0.870361664 0.901520531 1 0.916546112 0.897504521 0.955244123 0.870361664 0.901520531 Columns 7 through 12 0.97088762 0.915627954 0.88453985 0.976460857 0.872153443 0.923043135 0.97088762 0.915627954 0.88453985 0.976460857 0.872153443 0.923043135 Columns 13 through 18 1.025809862 0.941502448 0.904196842 0.984450007 0.859495215 0.916077493 1.025809862 0.941502448 0.904196842 0.984450007 0.859495215 0.916077493 Columns 19 through 24 0.997053985 0.912399374 0.878544039 0.969741146 0.885025052 0.926594672 0.997053985 0.912399374 0.878544039 0.969741146 0.885025052 0.926594672 Column 25 0.968193226 0.968193226
CheckDifference = diff(Check);
disp(CheckDifference)
Columns 1 through 6 0 -1.11022302462516e-16 0 0 0 1.11022302462516e-16 Columns 7 through 12 -1.11022302462516e-16 0 0 0 0 -1.11022302462516e-16 Columns 13 through 18 0 2.22044604925031e-16 0 1.11022302462516e-16 -1.11022302462516e-16 0 Columns 19 through 24 0 1.11022302462516e-16 1.11022302462516e-16 -1.11022302462516e-16 0 1.11022302462516e-16 Column 25 0
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
M = 0.0522
I = 9
P = 1/f(I); % Period [s]
P_hrs = P/3600 % Period [hrs]
P_hrs = 25
% vq2 = interp1(EXP_1(7,:),v,xq,'spline');
% plot(x,v,'o',xq,vq2,':.');
% % xlim([0 5E-11])
Xsize = numel(X)
Xsize = 25
% tinterp = linspace(min(tdata), max(tdata), N);
Xinterp = interp1(tdata, X, tinterp)
Xinterp = 1×49
1.0000 0.9583 0.9165 0.9070 0.8975 0.9264 0.9552 0.9128 0.8704 0.8859 0.9015 0.9362 0.9709 0.9433 0.9156 0.9001 0.8845 0.9305 0.9765 0.9243 0.8722 0.8976 0.9230 0.9744 1.0258 0.9837 0.9415 0.9228 0.9042 0.9443
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
format longG
Check = [X; Xinterp(1:2:end)];
disp(Check)
Columns 1 through 6 1 0.916546112 0.897504521 0.955244123 0.870361664 0.901520531 1 0.916546112 0.897504521 0.955244123 0.870361664 0.901520531 Columns 7 through 12 0.97088762 0.915627954 0.88453985 0.976460857 0.872153443 0.923043135 0.97088762 0.915627954 0.88453985 0.976460857 0.872153443 0.923043135 Columns 13 through 18 1.025809862 0.941502448 0.904196842 0.984450007 0.859495215 0.916077493 1.025809862 0.941502448 0.904196842 0.984450007 0.859495215 0.916077493 Columns 19 through 24 0.997053985 0.912399374 0.878544039 0.969741146 0.885025052 0.926594672 0.997053985 0.912399374 0.878544039 0.969741146 0.885025052 0.926594672 Column 25 0.968193226 0.968193226
CheckDifference = diff(Check);
disp(CheckDifference)
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
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 个评论
显示 4更早的评论
Ali
Ali 2024-12-3，15:47
Hello Star Strider!
I have two questions:
  • Can I provide boundaries for the frequency output of the fft analysis so that it returns values within the upper - lower limits of the specified frequency?
  • Why do you think my fft peaks are wider (the peaks are not as sharp as other frequency domain plots) than the general frequency domain plots in other references? Is it because of the low sampling rate I have?
Best,
Star Strider
Star Strider 2024-12-3，16:42
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
.

请先登录，再进行评论。

更多回答（2 个）

Sahas
Sahas 2024-11-1
Hi @Ali,
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:
I hope this is beneficial!
Paul
Paul 2024-11-2
编辑：Paul 2024-11-4
Ran in:
Hi Ali,
The doc page for interpft states:
"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.

类别

在 Help Center 和 File Exchange 中查找有关 Transforms 的更多信息

标签

产品


版本

R2024b

Community Treasure Hunt

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

Start Hunting!

Translated by