medfreq with delta-F threshold ?

I am using the medfreq function to extract fast freq changes (few ms) in a sinewave signal. The function seems to work pretty well and better than "tfridge" (less artifacts), "sst" (faster) and instfreq (less artifacts).

however, when the source causes the signal to change in amplitude and no freq changes are present, it is when my problem starts. If one would imagine a spectrogram rapresentation of this: the signal would be a constant line, with some fast peaks occurring here and there (chirps). Using the medfreq function these are nicely detected.

Changes in amplitude are caused by movements in the signal source and cause the spectrogram trace to have gaps. In absence of freq modulations (chirps), the medfreq function detects a lot of noise - in correspondence of those gaps.

This suggests me that the medfreq function works "point to point", otherwise the gaps would be detected as "sinks" in the medfreq trace (but I could be wrong). Is there a way to impose a fixed range of freq in which the median is searched ? or somehow make it so that gaps in the main freq component are not creating artifacts which would be confused for chirps ?

[Greg Dionne]

Would you be willing to post some sample data?

[LO]

Sure,

the first is a trace from a signal of relatively constant intensity in which freq modulations are present (and detected by medfreq, freq ranges from signal to +500Hz ca). Note: do not pay attention to the negative values of the trace on the side: it is the result of some postprocessing to average the peak values (is the result of subtraction of another medfreq trace from the one corresponding to the signal shown in the spectrogram...but this is not relevant here).

the other below is a signal of weaker and less constant amplitude (but same freq) with no chirps.the issue is more in this second figure: those peaks detected in the medfreq trace correspond to the gap visible in the spectrogram trace on the left

a similar issue would occur using other functions to determine the instant frequency of a signal: instfreq, tfridge, fsst and zerocrossing analysis. Comparing all those functions I found that - at least for my purposes - the medfreq function is more reliable and creates less artifacts... if not for this exception: when the signal has low power and has "gaps". In this case the errors are actually worse than those produced by these other functions.

[Greg Dionne]

In your first example it seems like there's just one fairly dominant frequency component. I would probably try instfreq (https://www.mathworks.com/help/signal/ref/instfreq.html). See https://www.mathworks.com/help/signal/ug/hilbert-transform-and-instantaneous-frequency.html for an example (assuming you have a relatively recent copy of the Signal Processing Toolbox).

Your second example seems to have other faint components which may (or may not) interfere with obtaining the frequency. If I couldn't get reliable results on that, then I'd probably try to downconvert to baseband and try to reconstruct the signal from there.

If you attach your data, I can give it a try.

[LO]

thanks Greg

see if you can open these mat files. they should be sinewave signals, half of those depicted above. In the second file there are some gaps. They both last 30 seconds . The sampling freq is 20kHz

Answer 1

[Greg Dionne]

From what I can tell your sinusoid is around ~895 Hz and has fairly clean second and third harmonics. So I took that as the starting point. The approach is to bandpass about each of the harmonics, mix down, filter and reconstruct the instantaneous frequency; if you divide each of the harmonics by its order, the overlaid results have reasonable agreement.

Anyways, this should hopefully get you started.

load seg1
Fs = 20e3;
%  carrier
Fc = 895;

% choose 100 Hz (one-sided) bandwidth about carrier
bw = 100;

% attempt FM demodulation about carrier
[Finst1, Tinst] = instfreq_bb(seg1, Fs, Fc, bw);
[Finst2, Tinst] = instfreq_bb(seg1, Fs, 2*Fc, bw);
[Finst3, Tinst] = instfreq_bb(seg1, Fs, 3*Fc, bw);

% superimpose 
plot(Tinst,[Finst1(:) Finst2(:)/2 Finst3(:)/3 (Finst1(:)+Finst2(:)/2+Finst3(:)/3)/3])
ylabel('Freq');
xlabel('Time');
legend('1st','2nd / 2','3rd / 3','average')

function [Finst, Tinst] = instfreq_bb(xx, Fs, Fc, bw)

%spectrogram(xx,kaiser(1024,10),1000,1024,20e3,'yaxis','power')

% pre-filter about carrier
x = bandpass(xx, Fc+[-bw bw], Fs, 'Steepness',.5);
%spectrogram(x,kaiser(1024,10),1000,1024,20e3,'yaxis','power')

t = (0:numel(x)-1)./Fs;

% mix down
z = complex(x .* cos(-2*pi*Fc*t), ...
            x .* sin(-2*pi*Fc*t));

% filter
bb = lowpass(z,bw,Fs,'Steepness',0.99);
%spectrogram(bb,kaiser(1024,10),1000,1024,20e3,'yaxis','power','centered')

% fetch instantaneous frequency from angular component.
Finst = angle(bb(2:end).*conj(bb(1:end-1))).*Fs/(2*pi)+Fc;

% fetch weighted time.
Tinst = t(1:end-1)+t(2)/2;
end

[LO]

Thanks Greg

I tried your script. could you please elaborate a bit more ? What would "mix down" mean ? I thought you were basically trying to reconstruct the signal from its basic components...by creating a new one which has the same harmonic frequencies. However plotting the data seem to be as fragmented as the inst freq trace. Maybe I missed something, sorry

[Greg Dionne]

Sure.

As I read your question it seemed like you were hoping to use medfreq to determine instantaneous frequency as you move along each column of a spectrogram, and you suspected changes in amplitude to be affecting your result. So the above code is an attempt to determine the instantaneous frequency of each component individually: if the overall signal is moving in frequency, then we should see all the harmonics of that signal move up or down in unison; when they are not, it may be due to some other form of noise.

Its possible I misinterpreted your setup though. How fast are the chirps you are expecting to see (a few ms, or tens of ms?) If there are unwanted amplitude components that are fast moving at a similar timescale, it may of course be difficult to separate out the two.

I also missed your question about selecting a range with medfreq. There is an optional FREQRANGE input which you can add at the end of your input list - that will restrict the range of frequencies searched. An example of that is below:

nSamp = 1024;
Fs = 1024e3;
t = (0:nSamp-1)'/Fs;
x = sin(2*pi*t*100.123e3);
 
[Pxx, F] = periodogram(x,kaiser(nSamp,38),[],Fs);
medfreq(Pxx,F,[50 150]*1e3)

You'll see an output like:

The shaded area is ingored for the purposes of computation.

Feel free to share the code you are using if you want me to take a look at it.

Hope this helps,

-Greg

[LO]

thanks Greg, I have tried to change the type of analysis window (from hanning to blackmanharris). I noticed this solved the issue (it has a stronger filtering for whatever is in the surrounding of the carrier fr).

By applying at the same time a threshold on the spectrogram it is possible to visualize a cleaner signal and calculate medfreq with better results.

However still I have some issues as I am not really able to select only those points in which the power is high enough.

in the formula for medfreq calculation one can set ranges for power detection based on freq values and time values. Perhaps if one could also impose a power threshold there.. this could solve the problem (the problem now is still the same: sometimes - in correspondence of power losses in the carrier fr - some peaks are detected. this might have something to do with the 3D nature of the spectrogram and to whatever medfreq detects in those spots... at least this is what I suspect. somehow maybe the solution would be to get some 2D plane out of that 3D mess).... am I right ?

[Greg Dionne]

If you are finding that changing the window has a favorable effect, you could try a parameterized window. Kaiser is often good (start with a low beta (say, 0.4) then gradually increase - but I think in your case, keep it under 10.

Eventualy though you will reach the limitations of what you can expect a median frequency method to do. It is simply looking for where the power level is split equally between the low and high frequency range. If there isn't enough discernable signal power, there's not much one can do...(the picture you posted looked like the power had completely disappeared into the noise.) I'm not sure even a reassigned spectrogram would be able to handle that case.

Perhaps other harmonics could contain enough power (or all of them together) to get a lock on your frequency? Maybe you could try computing the frequency of each harmonic separately, dividing by the harmonic number, and take weighted average by power level?

Hope this helps.

-Greg

[LO]

Thanks Greg

I forgot to mention I have tried to apply your suggestion: no, considering the harmonics does not make it any better. I agree with you: there is not much to do in case of noise.

however what I want to achieve is not measauring the signal in the noisy spots. Is the opposite (exclude it). this could be achieved a little bit by changing the window and by thresholding the spectrogram (without reassignment, as it would eventually increase the noise for medfreq). I was wondering whether there is another way to apply a threshold directly on medfreq (or whether a similar function exists).

well thanks for your feedback so far ! it was helpful