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

# FFT interpolation Matlab

Started by ●March 6, 2019

Reply by ●March 7, 20192019-03-07

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 >

Reply by ●March 7, 20192019-03-07

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

Reply by ●March 7, 20192019-03-07

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

Reply by ●March 7, 20192019-03-07

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