trigonometric interpolation of a sampled signal

Question

Given N sampled points, using the FFT we can get the Fourier transform of those N points $X_k$. With N/2 the Nyquist frequency and $X_0$ the DC value. Using the inverse we can then get back the original function we just measured. However if we would like more points then just the N we have measured but instead we would like M, how can u use the inverse FFT to find the trigonometric interpolation? We can assume the N is even and that M>N. And wat if we would drop values out of $X_k$, how would you find a trigonometric interpolation of the original signal.

Answer 1

What you want to do is treat the FFT of the signal as if the signal were sampled at a higher rate than it really was. You do this by simply padding the FFT of $x$ with zeros prior to taking the IFFT.

To understand this, consider some input signal $x[n]$ sampled at some frequency $f_s$. Since $x[n]$ is sampled at $f_s$, it has no frequency content less than $-f_s/2$ or greater than $+f_s/2$. This is true by definition, regardless of whether there was aliasing in the original sampling of $x$ or not. Also, we know that the FFT of $x[n]$ extends from $-f_s/2$ to $+f_s/2$. Since $x[n]$ has no frequency content outside of this range, what would its FFT look like if it were sampled at a higher rate? The answer is that the frequency bins of the FFT that were outside this range would have a value of zero. Thus, if we take the FFT of $x[n]$ and add zeros to it, it will look like the FFT of $x[n]$ but sampled at a higher rate. Hopefully, that makes sense.

One thing to be careful of is that you put the zeros in the right location. Most software does not shift the FFT. You want to add zeros at the locations that correspond to the higher frequency locations, which will be in the middle of an unshifted FFT. Hence, people usually swap the FFT, add zeros the left and right sides, unswap the new zeropadded array, and then perform the IFFT.

EDIT

Here's an example. Keep in mind that I do not have access to Matlab at the moment, so there might be some typos that I miss. I'm also using round numbers so that I can do some calculations in my head. In general, you'd want to write a function to figure out the amount of zeropadding required to achieve a desired sample rate and to do the actual zeropadding.

% First let us create a test signal. Quite arbitrarily I
% will us a 137(Hz) cosine sampled at 1(kHz).
fs = 1e3;                   % Sampling frequency (Hz).
t = 0 : 1/fs : 1-1/fs;      % 1000 time samples (s).
x = cos( 2*pi * 137 * t );  % Signal of interest.
N = length( x );            % Signal length.

% Now we take the FFT of x and zeropad it to make
% it look like it was sampled at 2(kHz). A 2(kHz) 
% sampling frequency would result in 2000 time
% samples and frequency bins rather than 1000.
% Thus, we want to add 1000 zeros to the FFT of x.
X = fftshift( fft( x ) );
Xz = [ zeros(1,500), X, zeros(1,500) ];
xup = ifft( ifftshift( Xz ) );

% Let's plot the two signals for comparison. They
% should be the same signal, just sampled at different
% rates.
fup = 2*fs;
tup = 0 : 1/fup : 1-1/fup;

figure;
plot( t, x, 'r' );
hold on;
plot( tup, xup, 'k' );

There may be a small residual imaginary component to xup that matlab complains about when you try to plot it. Just add a call to real in that case. The imaginary component is essentially numerical noise in this case.

Answer 2

Ok so if i understand correctly,

if i have for example x(n) points, lets say 10. And I chose M to be equal to 100. And i have an N samples (Note i am working with a matrix $N$x$d$, so i have multiple sampled signals)

N=size(x,1);
d=size(x,2);
for i=1:d
    Xk(:,i)=fft(x(:,i))
    for j=1:N
        if K<j && j<=N-K
             Xk(j,i)=0;
        end
    end
end

Here depending on the value of K it will drop values from the fft of the sampled signal. So basically $K<j<N-K$ will be equal to zero. once we have done that we can use the following relation to add "zeroes".

And we also also know that $Y_k=conjugate(Y_{M-k}$) for $k=M/2+1,....,M-1$

So once we have all those $Y_k$ we can basically use MATLAB ifft and we get the interpolation. Correct?

EDIT

Ok, but what i want to do is given a random points, which I draw myself, how can i find the best approximation of a function that connects those points. So basically an interpolation. And i also want to see the different harmonics, that is why i am dropping values of the FFT.

function [y,Xk,Y] = poly(x,K,M)
N=size(x,1);
d=size(x,2);
for i=1:d
    Xk(:,i)=fft(x(:,i)) %calculate the FFT of the x(n) pionts
    for j=K+1:N
             Xk(j,i)=0; %we will make the value of the Xk equal to zero for indexes  
                        %greater then K. So K index will be the the index of the last                                                                                                                        
                        %Xk coefficient. Basically the one with the highest frequency.
                        %After that we will make them equal to zero. Also note the when
                        %K=N/2 we will not drop any values!                
    end
end

%Here we will find the Yk using the provided formulas nul=zeros((M/2+1)-(N/2+2),1) Y=[(M/N).Xk(1:N/2);(M/N)(Xk(N/2+1)/2);nul;] conj=fliplr(Y'); Y=[Y;conj.']; y=ifft(Y); plot(y);hold on;plot(time,x(:,i)) end Is it correct that the data of the FFT after the Nyquist freq. is obsolete? So basically i can drop the data of the FFT after N/2 and still recreate the original function.

The above algorithm is based on the following. Assume that we want to evaluate $y(t)$ in M points, then we will want to calculate $y(t_m)$ where $t_m=m/M with m=0,...,M-1$.

however we cannot find the IFFT of the above formula that is why we will use

were

since M>N and M is even. That way we can find $(y_m)^{M-1}_{m=0}$ by using the IFFT on the $(Y_k)^{M-1}_{m=0}$

EDIT2

Its ok i found the solution, what i wated to do is the following. Basically the Gibbs-phenomena. I can use the same algorithm for other functions.