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