FFT result dependent on number of interpolation points
28 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Ali
il 1 Nov 2024
Commentato: Star Strider
il 3 Dic 2024 alle 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.
0 Commenti
Risposta accettata
Star Strider
il 1 Nov 2024
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 Commenti
Star Strider
il 3 Dic 2024 alle 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
.
Più risposte (2)
Sahas
il 1 Nov 2024
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!
Paul
il 2 Nov 2024
Modificato: Paul
il 4 Nov 2024
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.
0 Commenti
Vedere anche
Categorie
Scopri di più su Transforms in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!