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