Instantaneous Frequency Measurement

Parth VakilFebruary 4, 200821 comments

I would like to talk about the oft used method of measuring the carrier frequency in the world of Signal Collection and Characterization world. It is an elegant technique because of its simplicity. But, of course, with simplicity, there come drawbacks (sometimes...especially with this one!).

In the world of Radar detection and characterization, one of the key characteristics of interest is the carrier frequency of the signal. If the radar is pulsed, you will have a very wide bandwidth, a short burst but great SNR. If on the other hand, it is a spread spectrum radar, you will have a signal that is present for a long time but the SNR is low and hence processing and extracting information becomes a challange. In the first case, we need a technique that is able to give us an answer on the center frequency with very little information available. In the second case, we need to have a technique that is resiliant enough to operate in degraded conditions. IFM is one such technique. I will talk about the drawbacks in just a little bit.

Let us look at the basic mathematical foundation real quick:

Let us start with a simple example of a simple sinusoid with a constant amplitude.

s(t) = Acos(ωt + Φ)

where A is the constant amplitude, ω is the angular frequency and Φ is the initial phase.

If we were to delay this signal by some τ and multiply it with the original signal we will have:

s(t)s(t+τ) = Acos(ωt + Φ)Acos(ω(t+τ)+Φ)

s(t)s(t+τ) = A2⁄2·[cos(2ωt+ωτ+2Φ) - cos(-ωτ)]

If we low pas filter the above delay multiplied signal, we are left with:

s1(t) = -A2⁄2·cos(ωτ)

Now, if we were to measure the signal s(t) and s(t+τ) we can have the experimental value s1(t). If we perform the angle calculation (inverse cosine), we will end up with

ω = [arccos(-2s1(t)/A2)-Φ + 2kπ]/τ

Thing to note in digital systems is that we would typically have an analytic signal for processing rather than a real signal. This makes the arccos into an arctan. Still an expensive operation on a FPGA! So in the digital case, the s[n] we would be looking at would be an imaginary signal. If we follow though with the math above given s[n] = Aexp[ωn+Φ], we would arrive at a similar formula for ω above except the arccos will be substituted by an arctan.

Thing to note as well is that in either case (original signal real or complex) there is this 2 pi ambiguity. This will make our signal alias over the 2-pi interval. Hence we need to find a way to resolve this ambiguity. This is easily done by the observation that if we pick multiple delays (τ) and have multiple delay-multiply values, we can find the exact frequency. This is demonstrated through the following figures. This was done in MATLAB and you might find the attached script useful.

In the first two plots, the SNR is infinite. The sampling rate used is 1000 and the center frequency of the signal is 200. I have used three different delays (delay = 2, 4, 8). The plots of these delays are shown. Also shown is the output frequency of the IFM.

Phase Output for f = 200 and tau = 2, 4, 8

In the above plot, my intention with the arrows at the intersections is to show the ambiguity as τ increases. In short, while the resolution increases with τ, so does the ambiguity.

The graph of the frequency is:

Frequency Output

Now, if we look at the phase and frequency output for the same signal with a SNR of 2dB over the full bandwidth, following are the results:

2dB SNR Phase Outputs

And frequency:

Freuqnency Output at 2dB SNR

Some interesting features arise when the frequency we are trying to measure is related to the angles pi/4, pi/2, pi etc. because ωτ for various τ can give us a number that is either close to or exactly pi, 2pi, 0 etc. The issue with a number close to 0 or 2pi is the fact that the measured phase will jump between 0 and 2pi with noise. This will also ruin our frequency measurement. This is illustrated in the two figures below:

Frequency is 240, SNR = 2dB

Frequency Output

In the above two figures, the SNR is still 2dB but the frequency is now 240. The angular frequency for f = 240 is (240*2pi/1000). For delay of 4, ωτ ≈ 2pi. With noise, this value of ωτ will jump rapidly between 0 and 2pi as noticed in subplot 2 above. Also, for delay of 8, ωτ ≈ 4pi. Again, this value will jump between 0 and 2pi as shown in subplot 3 above. If you look at the frequency output now, it is quite noisy. A way to alleviate this problem is to unwrap the phase measurements.

The major drawback of the IFM scheme is presence of in band interferers. If the interferer is stronger than the signal we are trying to measure, the IFM will "pull" the result towards the stronger signal. In the digital realm, this is mitigated by having smaller bandwidths of interest and using channelizers.

The MATLAB code that I used to generate is as below:

fs = 1000;
fc = 240
t = 0:1/fs:4095/fs;
s = awgn(floor(127*cos(2*pi*fc*t)),2,'measured');
s_imag = hilbert(s);
d1 = 2;
d2 = 4;
d3 = 8;
lpf = [-2.6498e-019 7.977e-005 0.00036582 0.00095245 0.0019625 0.0035319 0.0057835 0.0087932 0.012558 0.01697 0.021811 0.026759 0.031423 0.035391 0.038282 0.039806 0.039806 0.038282 0.035391 0.031423 0.026759 0.021811 0.01697 0.012558 0.0087932 0.0057835 0.0035319 0.0019625 0.00095245 0.00036582 7.977e-005 -2.6498e-019];
s_d1 = [s_imag(d1+1:end) zeros(1,d1)];
s_d2 = [s_imag(d2+1:end) zeros(1,d2)];
s_d3 = [s_imag(d3+1:end) zeros(1,d3) ];

delay1 = angle(filter(lpf,1,s_d1.*conj(s_imag)));
delay2 = angle(filter(lpf,1,s_d2.*conj(s_imag)));
delay3 = angle(filter(lpf,1,s_d3.*conj(s_imag)));

delay1(delay1<=0) = delay1(delay1<=0) + 2*pi;
delay2(delay2<=0) = delay2(delay2<=0) + 2*pi;
delay3(delay3<=0) = delay3(delay3<=0) + 2*pi;

k1 = zeros(size(delay1));
k2 = zeros(size(delay2));

k1(delay1 > pi) = 2;
k2(delay2 > pi) = 1;

k = k1+k2;

w = (delay3 + 2*k*pi)/d3;
f = fs*w/(2*pi);

x = linspace(0,pi,size(delay3,2));
y1 = mod(2*x, 2*pi);
y2 = mod(4*x, 2*pi);
y3 = mod(8*x, 2*pi);
subplot(3,1,1); plot(x,delay1); hold on; plot(x, y1,'r'); title('TAU = 2'); grid on; set(gca, 'Ylim', [0 2*pi]);
subplot(3,1,2); plot(x,delay2); hold on; plot(x, y2,'r'); title('TAU = 4'); grid on; set(gca, 'Ylim', [0 2*pi]);
subplot(3,1,3); plot(x,delay3); hold on; plot(x, y3,'r'); title('TAU = 8'); grid on; set(gca, 'Ylim', [0 2*pi]);

plot(f); ylabel('Frequency'); title('Frequency Output of IFM')

You can play around with different sampling frequencies, center frequencies and SNRs to see how well the scheme works. You can also change the signal so that instead of being a simple sinusoid its a FM, phase coded or something else. I think, given the speed with which you arrive at an answer and relative implementation simplicity on a FPGA makes this scheme attractive.

I would like to add that please feel free to ask questions. If the concept of resolving the ambiguity is not clear, I can write up another short blog regarding how that works out. Hopefully you will be able to decipher the details from the MATLAB script. But if not, let me know.


1. "Digital Techniques for Wideband Receivers" by James Tsui

[ - ]
Comment by yhuangrFebruary 21, 2008
There shouldn't be a Φ in s1(t) = -A2⁄2·cos(ωτ+Φ), right?
[ - ]
Comment by vakilpFebruary 21, 2008
Rick: Good catch. Yes, the phi gets cancelled and the low passed version does not have the phi in it. Thanks. I made the change.
[ - ]
Comment by cpshah99March 23, 2008
Hi Parth Can you please tell me how to measure snr from the signal curropted by noise and multipath in real time? To mitigate multipath I am using adaptive DFE. Now I want to implement soft in soft out demapper(demodulator) and to implement that I need snr. Can I get the snr from PSD of the received signal?
[ - ]
Comment by vakilpMarch 23, 2008
Hi cpshah, I am not sure of how your system is set up. But what you could potentially do is terminate your antenna so that all your A/D is recording is the noise in your reciever. Since your SNR is mostly based on this, you can average over certain amount of time to figure out the standard deviation of the noise in the channel and hence the noise power. In your PSD, after you are recieving the signal, you can measure the max values or an average of the first 5-10 or more values to get an average signal power. Also, I am not sure of your constraints when you say real time. This method would take milli seconds for measurement. Your noise power ought not to be changing much over time unless your temperature is going to keep changing. Hope this helps. Parth
[ - ]
Comment by cpshah99March 23, 2008
HI Parth Thanks for reply. What I have done is that I designed coded signal in matlab and gave to my advisor for sea trial. Now what they did is that they transmitted the signal where transmitter was placed under the sea. And then they recorded the signal at receiver placed at 800m distance. Now what i have is recorded signal in .wav format. And to mitigate channel effects, I am using adaptive DFE equaliser. And my system is working with inbuilt demdulator of matlab but it gives hard o/p. I just want to implement soft demapper but the equation requires knowledge of noise variance. So I am stuck at this point. I tried to use welch's method but I am getting power values in -ve. Thanks for your help. Chintan
[ - ]
Comment by rashdkMay 9, 2008
hi parth , thanks for info content that has resolved some of my doubts..I tried ur matlab IFM code with my own test sample (generated in sound forge with noise superimposed over it). using wavread to read my sample wav file to determine its fundamental freq , do i need the "AWGN step of ur code ".. can u suggest changes in ur code to retrieve fundamental freq from my noise embedded sample with delays corresponding to reception of my sample from spatially dispersed mics. thanks n regards dk
[ - ]
Comment by vakilpMay 9, 2008
Hi dk, If the signal you are providing to the IFM code already has noise super imposed then you do not need the AWGN step. That simply adds noise to the signal in MATLAB. This code should successfully find the fundamental frequency that you are looking for. I am not sure what you mean by "spatially dispersed mics". My assumption is that you have multiple transducers and you are working with pressure waves. If you have multiple such recievers, you can use the same exact script except space your transducers apart at lambda/2, lambda and 2*lambda. You can then follow the math above to resolve the ambiguities that you find in lambda and 2*lambda spacings. A lot of research has been done on optimum spacing of antenna for angle of arrival measurement. All of what I have come across has focused on electromagnetic waves. With acoustic waves, you might have greater problems with multipath depending on your environment. For that, you can have an additional channel that does amplitude measurement to detect presence of the signal and then perform angle measurement. Hope this helps. Parth
[ - ]
Comment by rashdkMay 10, 2008
hi parth , thanks for ur guidance..i have been trying to do as u mentioned with the same matlab compilation error . The different lines of code used by me for my sample are : fs = 192000; s=wavread('E:\SAMPLES\test4.wav'); s_imag = hilbert(s); (rest follows ur code) the error in compilation is : "Error using ==> horzcat CAT arguments dimensions are not consistent. Error in ==> ifm1 at 9 s_d1 = [s_imag(d1+1:end) zeros(1,d1)];" also , i tried to figure out ur lpf values to no avail. I must admit i m not an expert in matlab , being a comp Sc student trying to solve acoustic detection prob by computational means ..the input reqd , however , is the fundamental freq where i m stuck for the moment. Apologies for consuming too much space here but a way ahead to resolve the above prob would do a world of good to my pace of work. thanks n regards dk
[ - ]
Comment by vakilpMay 10, 2008
dk, There are couple problems: 1. When you import the the signal s using wavread, maybe the vector is a M by 1 vector. If you perform the function size(s) after you read it in, you will find that it returns M by 1. My code is assuming a 1 by M vector. What you can do to resolve the horzcat issue is simply say s = s' i.e. s is equal to s transpose. This should solve your problem with the error. 2. I am not sure what you mean when you say "tried to figure out ur lpf values". Maybe you are trying to figure out the amplitude response of the filter? If you perform the function fvtool(lpf) in MATLAB, it will give you an amplitude vs. normalized frequency graph. You can also check out its phase response. Since it is a LPF, it will work with your sampling rate as well. Hope this helps. Parth
[ - ]
Comment by rashdkMay 11, 2008
hi parth thanks once again for ur input as the problem is resolved as far as execution is concerned. From Figure 2 that plots freq , how do i infer the fundamental freq from the rather cluttered plot in my case.. If I may digress for another opinion from u : As the actual code is in c++ , as of now , computation time for a given reference area (100m by 100m) is large to generate a correlation pattern , iteratively for each point. I strongly suspect weakness in my correlation algorithm in C++ code. Can u suggest any good correlation ( cross corr) algorithm that is easy on resources to compute. Any reading material that gives a threadbare details of such algorithm so that i can incorporate it in my code using mathematical knowledge. The detection part is achieved but i m stuck in visual representaion whereby i generate a video for an acoustic correlation (between transducers) in a given area for a given sound source (detected sound source). thanks n regards dk
[ - ]
Comment by vakilpMay 11, 2008
dk, IFM, as noted in the blog text, is not a great way to get the frequency in presence of other interferes. It sounds to me that that is the case with your simulated data. Is this correct? If so, the IFM will gravitate towards that strongest signal in your processing band width. One way to counter this is to build a channelizer. Cut your bandwidth down significantly and implement your detect and IFM scheme on each of the channels. Also, you detection algorithm could potentially help you out if you perform IFM after detection making sure only signals of your interest are being detected and you are able to filter the other garbage out. But the next step would be for you to build a channelizer. Cross correlation is going to slow going unless you reduce the number of points you have. You can break down your cross correlation into smaller squares that over lap. Even still, to generate the whole picture will take significant amount of time. You could look at FFT based techniques to generate cross correlation.
[ - ]
Comment by rashdkMay 12, 2008
hi parth.. i would keep the suggested point in mind to optimise detection. As of now ,i have envisaged two scenario of application , one being an open field (unmanned territory ) to detect intrusions and the other being underwater (which is expected to throw lot of DSP challenges )..As advised , i m studying FFT techniques to achieve cross corr.. As regards ur code , i just wanted to confirm the x-coordinate of ur IFM plot and can a specific variable be designated that computes freq of interest , in addition to visual inference from graph.. i will seek further guidance wrt DSP as n when i encounter some problems. Honestly , ur guidance has been a ready reckoner for DSP novice like me. thanks n regards dk
[ - ]
Comment by rashdkMay 28, 2008
hi parth At the very outset , my sincere gratitude for all ur inputs that helped me in getting my detection algorithm work with good accuracy. Having achieved an accurate fix , I intend representing the gamut of aocustic data in visual form to enable a console operator appreciate sound concentration from the display. As of now , I can think of waterfall display with freq , amplitude and time (distance from transducers ). I have the received wavefroms from my array of sensors. How do u suggest i should go about it or do u recommend alternate approach. thanks n regards dk
[ - ]
Comment by nitind2291April 7, 2009
Hello Parth, Although I am not an electronic guy but I do certain frequency measurements (10MHz) of my resonating crystal as a function of time, but unable to get a stable frequency, there is always great noise associated with the signals (+/- 40Hz). Although I make sure there is no disturbance form surroundings, still that does not work and noise is still there. Since for me even +/- 2-3 Hz shift is not tolerable, I do not know how to correct this problem. I am using an impedance analyzer to monitor resonant frequency of the quartz crystal. Please help. Thanks Nitin
[ - ]
Comment by farashiJune 23, 2009
hi i wonder if instantaneous frequncy can have negative value? i think because some derivative on phase will be used to calculate it so negative freuency mathematically is true. by the way i have some problem with the physical meaning of negative frequency!
[ - ]
Comment by ksk_dneFebruary 21, 2010
very informative
[ - ]
Comment by BurhemJuly 14, 2010
Hello, Some of the issues you raised are adressed in the following 2 papers, and also the book that provides a more recent update. Hope this helps. Regards. B. Boashash, "Estimating and Interpreting the Instantaneous Frequency of a Signal-Part I: Fundamentals", Proceedings of the IEEE, Vol. 80, No. 4, pp. 519-538, April 1992. B. Boashash, "Estimating and Interpreting the Instantaneous Frequency of a Signal-Part II: Algorithms", Proceedings of the IEEE, Vol. 80, No. 4, pp. 539-569, April 1992. B. Boashash, editor, “Time-Frequency Signal Analysis and Processing – A Comprehensive Reference”, Elsevier Science, Oxford, 2003, ISBN 0-08-044335-4.
[ - ]
Comment by juttakerenMarch 30, 2011
Hi Parth, is that measurement also can be applied for continuous wavelet transform method? I mean Hilbert transform replace with wavelet transform literally. Now I am working on wavelet transform and have difficulties to measure intantaneous frequency in Matlab. For information, I have done it with wavelet ridge method but it seem did not quite right. I need an advice. Thank you. Regards, Jutta
[ - ]
Comment by veena hegdeMay 12, 2013
It is told that If you can measure the Instantaneous frequency and pass it through wavelet filter tuned to fundamental frequency every time it is possible to extract system frequency variation .I need to extract heart's funadmental frequency
[ - ]
Comment by sandy3129August 11, 2016
has any one done the ifm reciever using xilinx system generator block sets rather than matlab , im not able to work out the phase lut block???
[ - ]
Comment by Atif_117May 3, 2018

Hi Parth i was following the book you mentioned and the research paper "FPGA-Based 1.2 GHz Bandwidth Digital Instantaneous Frequency Measurement Receiver" by James Helton, Chien-In Henry Chen, David M. Lin and James B. Y. Tsui, DOI 10.1109/ISQED.2008.162. This paper shows 2MHz RMS error at  max but according to your code at mutiples of pi, pi/2, pi/4 we get very large errors due to phase jumps. How is that possible to get 2MHz RMS Error performance in presence of these large errors.

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Please login (on the right) if you already have an account on this platform.

Otherwise, please use this form to register (free) an join one of the largest online community for Electrical/Embedded/DSP/FPGA/ML engineers: