FFT result dependent on number of interpolation points

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.

Answer 1

[Star Strider]

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.

.

[Ali]

I see, thank you for your clarification.

However, for different values of N (interpolarion points) I find the peak period (P_hrs) as half of it (N/2). Do you have any idea why fft result might depend on the number of interpolation points?

[Star Strider]

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]

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

[Star Strider]

My pleasure!

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

.

[Ali]

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]

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

.

Answer 2

[Sahas]

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:

• https://www.mathworks.com/help/matlab/ref/interpft.html
• https://www.mathworks.com/help/matlab/ref/double.interp1.html

I hope this is beneficial!

[Ali]

Thank you Sahas!

Answer 3

[Paul]

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.