10 views (last 30 days)

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
on 14 May 2019

Edited: Greg Dionne
on 14 May 2019

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

Greg Dionne
on 28 May 2019

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

Opportunities for recent engineering grads.

Apply TodayFind the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!
## 4 Comments

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/461294-medfreq-with-delta-f-threshold#comment_703988

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/461294-medfreq-with-delta-f-threshold#comment_703988

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/461294-medfreq-with-delta-f-threshold#comment_704173

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/461294-medfreq-with-delta-f-threshold#comment_704173

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/461294-medfreq-with-delta-f-threshold#comment_704456

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/461294-medfreq-with-delta-f-threshold#comment_704456

## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/461294-medfreq-with-delta-f-threshold#comment_704979

⋮## Direct link to this comment

https://www.mathworks.com/matlabcentral/answers/461294-medfreq-with-delta-f-threshold#comment_704979

Sign in to comment.