DSPRelated.com
Forums

Instantaneous frequency recovery using hilbert transform

Started by Yan.L October 3, 2004
Hello,

I am going to recover the instantaneous frequency f(t) of a chirp
signal y(t) using hilbert transform. Here is my MATLAB code:

Fs = 1000;		% Sampling frequency
T = 2;			% Signal time duration is 2 seconds.
t=0:1/Fs:T;                 
f0 = 50;                % Initial frequency is 50hz
f1 = 100;               % f(t) changes from 50hz to 100Hz in 2
seconds.

freq = f0 + (f1-f0)/T*t;

y = sin(2*pi*freq.*t);  % Chirp signal y(t)

figure(1); plot(t, freq); 

% Spectrum Plot of y(t)
N_FFT = 1024;           % Number of FFT points
y_FFT = fft(y,N_FFT);
f = round(Fs*(0:N_FFT/2)/N_FFT);           % Frequency scaling 
figure(2); plot(f, 20*log10(abs(y_FFT(1:N_FFT/2+1))))

% Hilbert tranform for instantaneous frequency recovery
output = real(y);
output_h = hilbert(output);
theta = angle(output_h);
ins_freq = angle(conv(ones(1,5),output_h(2:length(output_h)).*conj(output_h(1:length(output_h)-1))));
figure; plot((1:1:length(ins_freq))/Fs, ins_freq*Fs/2/pi)


The result I got is different from my initial setting. f(t) should
linearly change from 50hz to 100hz in 2seconds but in my result the
instantaneous frequency (ins_freq) changes (also linearly) from 50hz
to 150hz, i.e. the bandwidth is two times of my initial setting. Could
you please point out my error? Thanks a lot!
"Yan.L" <waters@starhub.net.sg> wrote in message
news:df850db4.0410031551.1f99b1ce@posting.google.com...

> Could > you please point out my error? Thanks a lot!
Your error is to assume that phase = 2*pi*frequency*time; the actual circumstance is that 2*pi*frequency = time derivative of phase. If f(t) = f0+(f1-f0)/(t1-t0)*(t-t0), then phase(t) = 2*pi*(f0+(f1-f0)/(t1-t0)*(t-t0)/2).*(t-t0)+phase0 -- write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, & 6.0134700243160014d-154/),(/'x'/)); end