FFT interpolation Matlab

Hi,
I am trying to use interpolation to detect an FFT peak between bins however my interpolation results in the same peak as the FFT.  I am inputting a 50Hz sine wave.. Sampling freq is 6500hz and I have 2048 sample points.  There fore my bin freq is 3.17hz which multipled out gives me a 50.8hz nearest bin frequency which is what my code is giving me.  I am trying to see where my interpolation code is going wrong and not giving me a peak at 50hz
Thanks in advance!

//Code:

Fs =6500;            % Sampling frequency                    
Y=SINWAVTest;
T = 1/Fs;             % Sampling period       
L = length(Y);             % Length of signal
t = (0:L-1)*T;        % Time vector
t2 = flipud(rot90(t));
y2 = Y; %+ 2*randn(size(t2));
y2 = y2;% - mean(y2);
y3 = fft((y2))
P2 = abs(y3/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
%%Interpolated 
T_mean = (t(L)-t(1))/(L-1);
t_u = t(1):T_mean:t(L);
yy = interp1(t,Y,t_u,'spline');
y4 = fft(yy)
Psd_i=abs(fft(yy)).^2/(L*T_mean)
fr=[-ceil((L-1)/2):floor((L-1)/2)]/(L*T_mean);
Psd_a=abs(nedft(Y,t2,fr)).^2/(L*T_mean);
P1I = Psd_i(1:L/2+1);
P1I(2:end-1) = 2*P1I(2:end-1);
%Plot time-domain signal
subplot(4,1,1);
plot(t, Y);
ylabel('Amplitude'); xlabel('Time (secs)');
axis tight;
title('Shaker Input Signal');
subplot(4,1,2);
plot(1000*t(1:1000),y2(1:1000))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')
subplot(4,1,3);
plot(f,P1)
 
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')
%ylim([0,2000]);
subplot(4,1,4);
plot(f,P1I,'r') 
title('Interpolated FFT')
xlabel('f (Hz)')
ylabel('|P1(f)|')
%ylim([0,2000]);

//Code

Well... I apologize if I say something wrong, but I hope I do not.

So... first, to put the terms right, what you are doing is DFT (discrete 
fourier transform), and FFT (fast fourier transform) is how you are 
doing it.

Second, when you do a DFT of a peak that doesn't exactly coincide with 
one of your DFT bins, you don't only get a single peak in that bin, but 
also a number of smaller peaks in other bins. Which is the source of 
your error. To reduce the number of bins containing any of the peaks, 
you need to multiply your signal by the window function prior to the 
DFT. That should massively reduce the error, and let you obtain much 
more accurate results.

Third, there's a bunch of window functions to choose from. I suggest you 
to start with something as simple as hamming window:
https://www.mathworks.com/help/signal/ref/hamming.html
After that you might try one of the more complicated windows such as 
Chebushev or Kaiser which allow you to control the magnitude of 
sidelobes, and see if it works better for you.

Fourth, at this point you might wonder, why does the lack of a window 
function affect your results at all? Without a window function (or, to 
put it in more precise terms, with a rectangular window function), you 
get noticeable sidelobes across all your DFT bins, but that shouldn't, 
in principle, affect the main peak, so you still should be able to get 
the correct result. Well, the problem here is that your signal contains 
not one complex sinusoid but two -- one at 50 Hz and the other at -50 
Hz. With a rectangular window function the sidelobes from either peak 
contaminate the other peak which results in a measurable error while 
determining its position.

Fifth, I hope my input was helpful, but you would be much better off if 
you read the first three chapters of "Understanding Digital Signal 
Processing" by Richard Lyons, that is the major source of my insight on 
the problem.

Hope that helps!

Gene.


On 07.03.2019 2:57, Michael R wrote:
> Hi,
> I am trying to use interpolation to detect an FFT peak between bins however my interpolation results in the same peak as the FFT.  I am inputting a 50Hz sine wave.. Sampling freq is 6500hz and I have 2048 sample points.  There fore my bin freq is 3.17hz which multipled out gives me a 50.8hz nearest bin frequency which is what my code is giving me.  I am trying to see where my interpolation code is going wrong and not giving me a peak at 50hz
> Thanks in advance!
>
> //Code:
>
> Fs =6500;            % Sampling frequency
> Y=SINWAVTest;
> T = 1/Fs;             % Sampling period
> L = length(Y);             % Length of signal
> t = (0:L-1)*T;        % Time vector
> t2 = flipud(rot90(t));
> y2 = Y; %+ 2*randn(size(t2));
> y2 = y2;% - mean(y2);
> y3 = fft((y2))
> P2 = abs(y3/L);
> P1 = P2(1:L/2+1);
> P1(2:end-1) = 2*P1(2:end-1);
> f = Fs*(0:(L/2))/L;
> %%Interpolated
> T_mean = (t(L)-t(1))/(L-1);
> t_u = t(1):T_mean:t(L);
> yy = interp1(t,Y,t_u,'spline');
> y4 = fft(yy)
> Psd_i=abs(fft(yy)).^2/(L*T_mean)
> fr=[-ceil((L-1)/2):floor((L-1)/2)]/(L*T_mean);
> Psd_a=abs(nedft(Y,t2,fr)).^2/(L*T_mean);
> P1I = Psd_i(1:L/2+1);
> P1I(2:end-1) = 2*P1I(2:end-1);
> %Plot time-domain signal
> subplot(4,1,1);
> plot(t, Y);
> ylabel('Amplitude'); xlabel('Time (secs)');
> axis tight;
> title('Shaker Input Signal');
> subplot(4,1,2);
> plot(1000*t(1:1000),y2(1:1000))
> title('Signal Corrupted with Zero-Mean Random Noise')
> xlabel('t (milliseconds)')
> ylabel('X(t)')
> subplot(4,1,3);
> plot(f,P1)
>
> title('Single-Sided Amplitude Spectrum of X(t)')
> xlabel('f (Hz)')
> ylabel('|P1(f)|')
> %ylim([0,2000]);
> subplot(4,1,4);
> plot(f,P1I,'r')
> title('Interpolated FFT')
> xlabel('f (Hz)')
> ylabel('|P1(f)|')
> %ylim([0,2000]);
>
> //Code
>

Also, I haven't seen you doing any interpolation in the code, which is 
why you might end up getting an integer multiple of a frequency bin as 
the output.

In that regard, it's essential to note that interpolation is to be done 
in the frequency domain, which could be achieved by either zero-padding 
of the time-domain series, or any conventional interpolation algorithms 
in the frequency domain (such as any kind of 3-point frequency 
estimation algorithms, which have been discussed extensively in this 
group in the past).

Best regards,
Gene.


On 07.03.2019 8:42, Gene Filatov wrote:
> Well... I apologize if I say something wrong, but I hope I do not.
>
> So... first, to put the terms right, what you are doing is DFT (discrete
> fourier transform), and FFT (fast fourier transform) is how you are
> doing it.
>
> Second, when you do a DFT of a peak that doesn't exactly coincide with
> one of your DFT bins, you don't only get a single peak in that bin, but
> also a number of smaller peaks in other bins. Which is the source of
> your error. To reduce the number of bins containing any of the peaks,
> you need to multiply your signal by the window function prior to the
> DFT. That should massively reduce the error, and let you obtain much
> more accurate results.
>
> Third, there's a bunch of window functions to choose from. I suggest you
> to start with something as simple as hamming window:
> https://www.mathworks.com/help/signal/ref/hamming.html
> After that you might try one of the more complicated windows such as
> Chebushev or Kaiser which allow you to control the magnitude of
> sidelobes, and see if it works better for you.
>
> Fourth, at this point you might wonder, why does the lack of a window
> function affect your results at all? Without a window function (or, to
> put it in more precise terms, with a rectangular window function), you
> get noticeable sidelobes across all your DFT bins, but that shouldn't,
> in principle, affect the main peak, so you still should be able to get
> the correct result. Well, the problem here is that your signal contains
> not one complex sinusoid but two -- one at 50 Hz and the other at -50
> Hz. With a rectangular window function the sidelobes from either peak
> contaminate the other peak which results in a measurable error while
> determining its position.
>
> Fifth, I hope my input was helpful, but you would be much better off if
> you read the first three chapters of "Understanding Digital Signal
> Processing" by Richard Lyons, that is the major source of my insight on
> the problem.
>
> Hope that helps!
>
> Gene.
>
>
> On 07.03.2019 2:57, Michael R wrote:
>> Hi,
>> I am trying to use interpolation to detect an FFT peak between bins
>> however my interpolation results in the same peak as the FFT.  I am
>> inputting a 50Hz sine wave.. Sampling freq is 6500hz and I have 2048
>> sample points.  There fore my bin freq is 3.17hz which multipled out
>> gives me a 50.8hz nearest bin frequency which is what my code is
>> giving me.  I am trying to see where my interpolation code is going
>> wrong and not giving me a peak at 50hz
>> Thanks in advance!
>>
>> //Code:
>>
>> Fs =6500;            % Sampling frequency
>> Y=SINWAVTest;
>> T = 1/Fs;             % Sampling period
>> L = length(Y);             % Length of signal
>> t = (0:L-1)*T;        % Time vector
>> t2 = flipud(rot90(t));
>> y2 = Y; %+ 2*randn(size(t2));
>> y2 = y2;% - mean(y2);
>> y3 = fft((y2))
>> P2 = abs(y3/L);
>> P1 = P2(1:L/2+1);
>> P1(2:end-1) = 2*P1(2:end-1);
>> f = Fs*(0:(L/2))/L;
>> %%Interpolated
>> T_mean = (t(L)-t(1))/(L-1);
>> t_u = t(1):T_mean:t(L);
>> yy = interp1(t,Y,t_u,'spline');
>> y4 = fft(yy)
>> Psd_i=abs(fft(yy)).^2/(L*T_mean)
>> fr=[-ceil((L-1)/2):floor((L-1)/2)]/(L*T_mean);
>> Psd_a=abs(nedft(Y,t2,fr)).^2/(L*T_mean);
>> P1I = Psd_i(1:L/2+1);
>> P1I(2:end-1) = 2*P1I(2:end-1);
>> %Plot time-domain signal
>> subplot(4,1,1);
>> plot(t, Y);
>> ylabel('Amplitude'); xlabel('Time (secs)');
>> axis tight;
>> title('Shaker Input Signal');
>> subplot(4,1,2);
>> plot(1000*t(1:1000),y2(1:1000))
>> title('Signal Corrupted with Zero-Mean Random Noise')
>> xlabel('t (milliseconds)')
>> ylabel('X(t)')
>> subplot(4,1,3);
>> plot(f,P1)
>>
>> title('Single-Sided Amplitude Spectrum of X(t)')
>> xlabel('f (Hz)')
>> ylabel('|P1(f)|')
>> %ylim([0,2000]);
>> subplot(4,1,4);
>> plot(f,P1I,'r')
>> title('Interpolated FFT')
>> xlabel('f (Hz)')
>> ylabel('|P1(f)|')
>> %ylim([0,2000]);
>>
>> //Code
>>
>

Thank you for the reply Gene,  I have ordered the Lyons book and it should be here in a couple days. Based on your direction I was able to figure out my code to work..

Thanks again!

Mike

On Thursday, March 7, 2019 at 3:59:03 AM UTC-8, Gene Filatov wrote:
> Also, I haven't seen you doing any interpolation in the code, which is 
> why you might end up getting an integer multiple of a frequency bin as 
> the output.
> 
> In that regard, it's essential to note that interpolation is to be done 
> in the frequency domain, which could be achieved by either zero-padding 
> of the time-domain series, or any conventional interpolation algorithms 
> in the frequency domain (such as any kind of 3-point frequency 
> estimation algorithms, which have been discussed extensively in this 
> group in the past).
> 
> Best regards,
> Gene.
> 
> 
> On 07.03.2019 8:42, Gene Filatov wrote:
> > Well... I apologize if I say something wrong, but I hope I do not.
> >
> > So... first, to put the terms right, what you are doing is DFT (discrete
> > fourier transform), and FFT (fast fourier transform) is how you are
> > doing it.
> >
> > Second, when you do a DFT of a peak that doesn't exactly coincide with
> > one of your DFT bins, you don't only get a single peak in that bin, but
> > also a number of smaller peaks in other bins. Which is the source of
> > your error. To reduce the number of bins containing any of the peaks,
> > you need to multiply your signal by the window function prior to the
> > DFT. That should massively reduce the error, and let you obtain much
> > more accurate results.
> >
> > Third, there's a bunch of window functions to choose from. I suggest you
> > to start with something as simple as hamming window:
> > https://www.mathworks.com/help/signal/ref/hamming.html
> > After that you might try one of the more complicated windows such as
> > Chebushev or Kaiser which allow you to control the magnitude of
> > sidelobes, and see if it works better for you.
> >
> > Fourth, at this point you might wonder, why does the lack of a window
> > function affect your results at all? Without a window function (or, to
> > put it in more precise terms, with a rectangular window function), you
> > get noticeable sidelobes across all your DFT bins, but that shouldn't,
> > in principle, affect the main peak, so you still should be able to get
> > the correct result. Well, the problem here is that your signal contains
> > not one complex sinusoid but two -- one at 50 Hz and the other at -50
> > Hz. With a rectangular window function the sidelobes from either peak
> > contaminate the other peak which results in a measurable error while
> > determining its position.
> >
> > Fifth, I hope my input was helpful, but you would be much better off if
> > you read the first three chapters of "Understanding Digital Signal
> > Processing" by Richard Lyons, that is the major source of my insight on
> > the problem.
> >
> > Hope that helps!
> >
> > Gene.
> >
> >
> > On 07.03.2019 2:57, Michael R wrote:
> >> Hi,
> >> I am trying to use interpolation to detect an FFT peak between bins
> >> however my interpolation results in the same peak as the FFT.  I am
> >> inputting a 50Hz sine wave.. Sampling freq is 6500hz and I have 2048
> >> sample points.  There fore my bin freq is 3.17hz which multipled out
> >> gives me a 50.8hz nearest bin frequency which is what my code is
> >> giving me.  I am trying to see where my interpolation code is going
> >> wrong and not giving me a peak at 50hz
> >> Thanks in advance!
> >>
> >> //Code:
> >>
> >> Fs =6500;            % Sampling frequency
> >> Y=SINWAVTest;
> >> T = 1/Fs;             % Sampling period
> >> L = length(Y);             % Length of signal
> >> t = (0:L-1)*T;        % Time vector
> >> t2 = flipud(rot90(t));
> >> y2 = Y; %+ 2*randn(size(t2));
> >> y2 = y2;% - mean(y2);
> >> y3 = fft((y2))
> >> P2 = abs(y3/L);
> >> P1 = P2(1:L/2+1);
> >> P1(2:end-1) = 2*P1(2:end-1);
> >> f = Fs*(0:(L/2))/L;
> >> %%Interpolated
> >> T_mean = (t(L)-t(1))/(L-1);
> >> t_u = t(1):T_mean:t(L);
> >> yy = interp1(t,Y,t_u,'spline');
> >> y4 = fft(yy)
> >> Psd_i=abs(fft(yy)).^2/(L*T_mean)
> >> fr=[-ceil((L-1)/2):floor((L-1)/2)]/(L*T_mean);
> >> Psd_a=abs(nedft(Y,t2,fr)).^2/(L*T_mean);
> >> P1I = Psd_i(1:L/2+1);
> >> P1I(2:end-1) = 2*P1I(2:end-1);
> >> %Plot time-domain signal
> >> subplot(4,1,1);
> >> plot(t, Y);
> >> ylabel('Amplitude'); xlabel('Time (secs)');
> >> axis tight;
> >> title('Shaker Input Signal');
> >> subplot(4,1,2);
> >> plot(1000*t(1:1000),y2(1:1000))
> >> title('Signal Corrupted with Zero-Mean Random Noise')
> >> xlabel('t (milliseconds)')
> >> ylabel('X(t)')
> >> subplot(4,1,3);
> >> plot(f,P1)
> >>
> >> title('Single-Sided Amplitude Spectrum of X(t)')
> >> xlabel('f (Hz)')
> >> ylabel('|P1(f)|')
> >> %ylim([0,2000]);
> >> subplot(4,1,4);
> >> plot(f,P1I,'r')
> >> title('Interpolated FFT')
> >> xlabel('f (Hz)')
> >> ylabel('|P1(f)|')
> >> %ylim([0,2000]);
> >>
> >> //Code
> >>
> >

Mike, you are welcome!

Gene.


On 07.03.2019 22:57, Michael R wrote:
> Thank you for the reply Gene,  I have ordered the Lyons book and it should be here in a couple days. Based on your direction I was able to figure out my code to work..
>
> Thanks again!
>
> Mike
>
> On Thursday, March 7, 2019 at 3:59:03 AM UTC-8, Gene Filatov wrote:
>> Also, I haven't seen you doing any interpolation in the code, which is
>> why you might end up getting an integer multiple of a frequency bin as
>> the output.
>>
>> In that regard, it's essential to note that interpolation is to be done
>> in the frequency domain, which could be achieved by either zero-padding
>> of the time-domain series, or any conventional interpolation algorithms
>> in the frequency domain (such as any kind of 3-point frequency
>> estimation algorithms, which have been discussed extensively in this
>> group in the past).
>>
>> Best regards,
>> Gene.
>>
>>
>> On 07.03.2019 8:42, Gene Filatov wrote:
>>> Well... I apologize if I say something wrong, but I hope I do not.
>>>
>>> So... first, to put the terms right, what you are doing is DFT (discrete
>>> fourier transform), and FFT (fast fourier transform) is how you are
>>> doing it.
>>>
>>> Second, when you do a DFT of a peak that doesn't exactly coincide with
>>> one of your DFT bins, you don't only get a single peak in that bin, but
>>> also a number of smaller peaks in other bins. Which is the source of
>>> your error. To reduce the number of bins containing any of the peaks,
>>> you need to multiply your signal by the window function prior to the
>>> DFT. That should massively reduce the error, and let you obtain much
>>> more accurate results.
>>>
>>> Third, there's a bunch of window functions to choose from. I suggest you
>>> to start with something as simple as hamming window:
>>> https://www.mathworks.com/help/signal/ref/hamming.html
>>> After that you might try one of the more complicated windows such as
>>> Chebushev or Kaiser which allow you to control the magnitude of
>>> sidelobes, and see if it works better for you.
>>>
>>> Fourth, at this point you might wonder, why does the lack of a window
>>> function affect your results at all? Without a window function (or, to
>>> put it in more precise terms, with a rectangular window function), you
>>> get noticeable sidelobes across all your DFT bins, but that shouldn't,
>>> in principle, affect the main peak, so you still should be able to get
>>> the correct result. Well, the problem here is that your signal contains
>>> not one complex sinusoid but two -- one at 50 Hz and the other at -50
>>> Hz. With a rectangular window function the sidelobes from either peak
>>> contaminate the other peak which results in a measurable error while
>>> determining its position.
>>>
>>> Fifth, I hope my input was helpful, but you would be much better off if
>>> you read the first three chapters of "Understanding Digital Signal
>>> Processing" by Richard Lyons, that is the major source of my insight on
>>> the problem.
>>>
>>> Hope that helps!
>>>
>>> Gene.
>>>
>>>
>>> On 07.03.2019 2:57, Michael R wrote:
>>>> Hi,
>>>> I am trying to use interpolation to detect an FFT peak between bins
>>>> however my interpolation results in the same peak as the FFT.  I am
>>>> inputting a 50Hz sine wave.. Sampling freq is 6500hz and I have 2048
>>>> sample points.  There fore my bin freq is 3.17hz which multipled out
>>>> gives me a 50.8hz nearest bin frequency which is what my code is
>>>> giving me.  I am trying to see where my interpolation code is going
>>>> wrong and not giving me a peak at 50hz
>>>> Thanks in advance!
>>>>
>>>> //Code:
>>>>
>>>> Fs =6500;            % Sampling frequency
>>>> Y=SINWAVTest;
>>>> T = 1/Fs;             % Sampling period
>>>> L = length(Y);             % Length of signal
>>>> t = (0:L-1)*T;        % Time vector
>>>> t2 = flipud(rot90(t));
>>>> y2 = Y; %+ 2*randn(size(t2));
>>>> y2 = y2;% - mean(y2);
>>>> y3 = fft((y2))
>>>> P2 = abs(y3/L);
>>>> P1 = P2(1:L/2+1);
>>>> P1(2:end-1) = 2*P1(2:end-1);
>>>> f = Fs*(0:(L/2))/L;
>>>> %%Interpolated
>>>> T_mean = (t(L)-t(1))/(L-1);
>>>> t_u = t(1):T_mean:t(L);
>>>> yy = interp1(t,Y,t_u,'spline');
>>>> y4 = fft(yy)
>>>> Psd_i=abs(fft(yy)).^2/(L*T_mean)
>>>> fr=[-ceil((L-1)/2):floor((L-1)/2)]/(L*T_mean);
>>>> Psd_a=abs(nedft(Y,t2,fr)).^2/(L*T_mean);
>>>> P1I = Psd_i(1:L/2+1);
>>>> P1I(2:end-1) = 2*P1I(2:end-1);
>>>> %Plot time-domain signal
>>>> subplot(4,1,1);
>>>> plot(t, Y);
>>>> ylabel('Amplitude'); xlabel('Time (secs)');
>>>> axis tight;
>>>> title('Shaker Input Signal');
>>>> subplot(4,1,2);
>>>> plot(1000*t(1:1000),y2(1:1000))
>>>> title('Signal Corrupted with Zero-Mean Random Noise')
>>>> xlabel('t (milliseconds)')
>>>> ylabel('X(t)')
>>>> subplot(4,1,3);
>>>> plot(f,P1)
>>>>
>>>> title('Single-Sided Amplitude Spectrum of X(t)')
>>>> xlabel('f (Hz)')
>>>> ylabel('|P1(f)|')
>>>> %ylim([0,2000]);
>>>> subplot(4,1,4);
>>>> plot(f,P1I,'r')
>>>> title('Interpolated FFT')
>>>> xlabel('f (Hz)')
>>>> ylabel('|P1(f)|')
>>>> %ylim([0,2000]);
>>>>
>>>> //Code
>>>>
>>>
>