MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

New to MATLAB?

# Thread Subject: Non stationary FFT problem - strange HELP!

Subject: Non stationary FFT problem - strange HELP!

From: Dean

### Dean (view profile)

Date: 12 Feb, 2010 13:12:04

Message: 1 of 3

I've a time domain signal defined as:

y = cos(2*pi* (f0-a*t) t);

where f0=10 and a=0.1, I took 4096 data samples
in the total period T=20 sec, which the sampling frequency is

f0 = 10;
N = 2^12;
t = linspace(0,20,N); NFFT = 2^nextpow2(N);
fs = 1/(t(2)-t(1));
a = 0.1;
ff = f0 - a*t;
y = cos(2*pi*ff.*t);
Y = fft(y,NFFT)/N; f = fs/2*linspace(0,1,NFFT/2+1);
subplot(221); plot(t,ff); ylim([5,12]); text(2,9.2,['\alpha=',num2str(a)]);
title('Beat frequency');
subplot(222); plot(t,y,'r-'); title('Time-Domain sig.')
subplot(212);
plot(f,20*log10(2*abs(Y(1:NFFT/2+1)))); xlim([5,15]); ylim([-50,2]);
% plot(f,2*abs(Y(1:NFFT/2+1))); % xlim([0,15]); ylim([-50,2]);
title('Power Spectrum')
vline(f0); grid on;

Wired thing is that the spectrum spam from 6Hz to 10Hz, and I know there was no frequency contents below 8Hz, f0-a*t=10-0.1*20=8.

Can you please check if there is any problem with my code, or could you please confirm if my code was actually correct, i.e. the FT is no good for this kind of scenarios.

Cheers!

Subject: Non stationary FFT problem - strange HELP!

Date: 12 Feb, 2010 14:55:26

Message: 2 of 3

Hi Dean,

I would like to remind you the following:

1. When you want to plot the non-negative side of the spectrum, your frequency values are (0:NFFT/2-1)*fs/NFFT.

2. Discrete fourier transform is closely related to Fourier series. Therefore, unless your signal is truly periodic [that is, the last sample of your signal is expected to be followed by the very first sample of your signal, e.g. sine of a single frequency and of full period length], you will see a spectrum that is not what you expect to see.

You can easily verify that your original signal is not periodic. So, it is not possible to see the spectrum you want.

3. Imagine that you are copying and pasting your signal before and after your original vector [This is what DFT assumes. You have circular convolution and stuff, remember?]. You can observe that there is an abrupt frequency change at the points of connection [8 to 10 Hz]. This will also wake up some of your spectral content.

Best.

"Dean " <wj.antispam@gmail.com> wrote in message <hl3k34$p3r$1@fred.mathworks.com>...
> I've a time domain signal defined as:
>
> y = cos(2*pi* (f0-a*t) t);
>
> where f0=10 and a=0.1, I took 4096 data samples
> in the total period T=20 sec, which the sampling frequency is
>
> f0 = 10;
> N = 2^12;
> t = linspace(0,20,N); NFFT = 2^nextpow2(N);
> fs = 1/(t(2)-t(1));
> a = 0.1;
> ff = f0 - a*t;
> y = cos(2*pi*ff.*t);
> Y = fft(y,NFFT)/N; f = fs/2*linspace(0,1,NFFT/2+1);
> subplot(221); plot(t,ff); ylim([5,12]); text(2,9.2,['\alpha=',num2str(a)]);
> title('Beat frequency');
> subplot(222); plot(t,y,'r-'); title('Time-Domain sig.')
> subplot(212);
> plot(f,20*log10(2*abs(Y(1:NFFT/2+1)))); xlim([5,15]); ylim([-50,2]);
> % plot(f,2*abs(Y(1:NFFT/2+1))); % xlim([0,15]); ylim([-50,2]);
> title('Power Spectrum')
> vline(f0); grid on;
>
> Wired thing is that the spectrum spam from 6Hz to 10Hz, and I know there was no frequency contents below 8Hz, f0-a*t=10-0.1*20=8.
>
> Can you please check if there is any problem with my code, or could you please confirm if my code was actually correct, i.e. the FT is no good for this kind of scenarios.
>
> Cheers!

Subject: Non stationary FFT problem - strange HELP!

From: Wayne King

### Wayne King (view profile)

Date: 12 Feb, 2010 19:28:22

Message: 3 of 3

"Dean " <wj.antispam@gmail.com> wrote in message <hl3k34$p3r$1@fred.mathworks.com>...
> I've a time domain signal defined as:
>
> y = cos(2*pi* (f0-a*t) t);
>
> where f0=10 and a=0.1, I took 4096 data samples
> in the total period T=20 sec, which the sampling frequency is
>
> f0 = 10;
> N = 2^12;
> t = linspace(0,20,N); NFFT = 2^nextpow2(N);
> fs = 1/(t(2)-t(1));
> a = 0.1;
> ff = f0 - a*t;
> y = cos(2*pi*ff.*t);
> Y = fft(y,NFFT)/N; f = fs/2*linspace(0,1,NFFT/2+1);
> subplot(221); plot(t,ff); ylim([5,12]); text(2,9.2,['\alpha=',num2str(a)]);
> title('Beat frequency');
> subplot(222); plot(t,y,'r-'); title('Time-Domain sig.')
> subplot(212);
> plot(f,20*log10(2*abs(Y(1:NFFT/2+1)))); xlim([5,15]); ylim([-50,2]);
> % plot(f,2*abs(Y(1:NFFT/2+1))); % xlim([0,15]); ylim([-50,2]);
> title('Power Spectrum')
> vline(f0); grid on;
>
> Wired thing is that the spectrum spam from 6Hz to 10Hz, and I know there was no frequency contents below 8Hz, f0-a*t=10-0.1*20=8.
>
> Can you please check if there is any problem with my code, or could you please confirm if my code was actually correct, i.e. the FT is no good for this kind of scenarios.
>
> Cheers!

Hi Dean, I think you may have calculated something wrong, because I get the minimum instantaneous frequency to be 6 Hz and the maximum to be 10 Hz which agrees with the plot.

Let's look at your signal code:

f0 = 10;
N = 2^12;
t = linspace(0,20,N);
Fs = 2^12/20;
a = 0.1;
ff = f0 - a*t;
y = cos(2*pi*ff.*t);

Writing your FM signal as y(t) = cos(theta(t)) with theta(t)=2*pi*10*t-0.1*2*pi*t.^2;
By the way, did you realize that you multiplying your a*t^2 term by 2*pi as well? At any rate the instantaneous frequency is 1/2*pi times the derivative of your theta(t), which I calculate to be 10-0.2*t so your instantaneous frequencies are:

instfreq = 10-0.2*t;
plot(t,instfreq)
min(instfreq)
max(instfreq)

Which agrees with the spectrum:

H = spectrum.periodogram;
plot(psd(H,y,'Fs',Fs)); axis([4 12 -40 0]);

Hope that helps,
Wayne