FFT result dependent on number of interpolation points

Question

Ali 2024-11-1，3:19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2163390-fft-result-dependent-on-number-of-interpolation-points

编辑： Paul 2024-11-4，13:16

FFT_96hrs.m

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Star Strider 2024-11-1，11:13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2163390-fft-result-dependent-on-number-of-interpolation-points#answer_1539875

在 MATLAB Online 中打开

FFT_96hrs.m

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>

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.

.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Star Strider 2024-11-1，14:33

在 MATLAB Online 中打开

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)

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)

Y =

Columns 1 through 9 0.0000 + 0.0000i -0.1076 + 0.1006i 0.3314 + 0.0254i 0.0266 + 0.0266i 0.5699 - 0.0364i -0.0934 + 0.0945i 0.0549 - 0.0297i 0.1269 - 0.0399i 1.0960 + 0.7061i Columns 10 through 18 -0.3055 - 0.1947i 0.0107 - 0.0551i 0.0226 - 0.0824i -0.1125 - 0.3007i -0.0062 - 0.0693i 0.0051 - 0.0396i 0.0089 - 0.0270i 0.0106 - 0.0199i 0.0115 - 0.0152i Columns 19 through 27 0.0120 - 0.0118i 0.0122 - 0.0092i 0.0124 - 0.0071i 0.0125 - 0.0053i 0.0125 - 0.0036i 0.0125 - 0.0021i 0.0126 - 0.0007i 0.0126 + 0.0007i 0.0125 + 0.0021i Columns 28 through 36 0.0125 + 0.0036i 0.0125 + 0.0053i 0.0124 + 0.0071i 0.0122 + 0.0092i 0.0120 + 0.0118i 0.0115 + 0.0152i 0.0106 + 0.0199i 0.0089 + 0.0270i 0.0051 + 0.0396i Columns 37 through 45 -0.0062 + 0.0693i -0.1125 + 0.3007i 0.0226 + 0.0824i 0.0107 + 0.0551i -0.3055 + 0.1947i 1.0960 - 0.7061i 0.1269 + 0.0399i 0.0549 + 0.0297i -0.0934 - 0.0945i Columns 46 through 49 0.5699 + 0.0364i 0.0266 - 0.0266i 0.3314 - 0.0254i -0.1076 - 0.1006i

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))

Fs2 = 1.3889e-04

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

M = 0.0522

I = 9

P = 1/f(I); % Period [s]

P_hrs = P/3600 % Period [hrs]

P_hrs = 12.5000

% 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 2024-11-2，2:53

Thank you Star Strider! Your answer helped me a lot.

Star Strider 2024-11-2，10:07

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

请先登录，再进行评论。

Answer 2

Sahas 2024-11-1，6:58

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2163390-fft-result-dependent-on-number-of-interpolation-points#answer_1539805

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!

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Ali 2024-11-2，2:52

Thank you Sahas!

请先登录，再进行评论。

Answer 3

Paul 2024-11-2，15:16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2163390-fft-result-dependent-on-number-of-interpolation-points#answer_1540100

编辑：Paul 2024-11-4，13:16

在 MATLAB Online 中打开

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

FFT result dependent on number of interpolation points

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（2 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

FFT result dependent on number of interpolation points

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（2 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论