Non stationary FFT problem - strange HELP!

[thedean515]

The time domain signal is defined as:\[
y\left(t\right)=A\cos\left(2\pi\left(f_{0}-\alpha t\right)t\right),\]
where A=1, \[f_{0}=10\] and \[\alpha=0.1\], I take 4096 data samples
in the total period T=20 sec, which the sampling frequency is
about 205 Hz.

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;



Am I right in thinking that there are shouldn't be any frequency contents below 8 Hz, but why the spectrum contradict the thinking, moreover, the spectrum is oddly symmetrical.

Would you please help me to find out where I have done wrong or can you confirm this is the correct numerical simulation to this problem.

Cheers!![/img]

 

[albbg]

No, you thinking is not correct. The argument of a trigonometric funcion represent a phase, not a frequency. Since the frequency is the first derivative of the phase with respect to the time, divided by 1/(2π), you can find the frequency applying:

freq = d/dt[2*π*(fo - α*t)*t]/(2*π) = [2*π*fo - 2*2*π*α*t]/(2*π)

then

freq = fo - 2*α*t

so, if fo = 10 Hz and α = 0.1 you will find:

freq(min) = 6 Hz and freq(max) = 10 Hz

Folding of the spectrum has to be studied if fo < 2*α*t;

Just a question: where you time signal is coming from (just curiosity) ?

Regards
Alberto

 

[thedean515]

Dear albbg,

Thank you so much for prompt reply and it was extremely useful!

I totally agree with you, now the spectrum makes more sense.

Is it possible to estimate the phase of the signal if only the discrete time-domain signal is known? I suppose I can start with the spectrum estimated, but the only method I can think of is to approximate the area under the power spectrum as:


\[\int_0^T y(t) dt\]

This doesn't look right, what do you think?

Many thanks!

 

[albbg]

The formula you propose to estimate the phase is just the average of y(t) multiplied by T, so it is not correct.
If the analytical form of your signal is that you shown, I think you could just recover the phase applying acos[y(t)] after normalization of your data stream in the range [-1, +1].
A constant phase offset will be added to your phase estimation since the phase will start from zero at the first sampling instant for which t0 = 0.
Furthermore your phase estimation will be wrapped in the range [0,2*Π]