How can I find and plot the change in frequency of a signal?

I have the following signal:

I'm looking for a way to find the change in the frequency and then plot it (ideally, I would like to find the equation that represents the change in frequency of the signal).

I've tried using the spectrogram function, along with various methods using ffts but I can't seem to figure out a way to do this at all. Everything seems to be dedicated to finding variations in amplitude or deconstructing the signal into its basic frequency components (as with a fourier transform). Is such an analysis possible, and if so how might I go about doing this?

As always, thank you in advance for your help.

-Chris

Answer 1

[Star Strider]

The easiest way I can think of is to find the zero-crossings of the signal and take the inverse to calculate the frequency.

Example:

x = linspace(0,10*pi,1E+5);                                             % Create Data
y = sin(1./(x+0.01));                                                   % Create Data

zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                    % Returns Approximate Zero-Crossing Indices Of Argument Vector

dy = zci(y);                                                            % Indices of Approximate Zero-Crossings
for k1 = 1:size(dy,1)-1
    b = [[1;1] [x(dy(k1)); x(dy(k1)+1)]]\[y(dy(k1)); y(dy(k1)+1)];      % Linear Fit Near Zero-Crossings
    x0(k1) = -b(1)/b(2);                                                % Interpolate ‘Exact’ Zero Crossing
    mb(:,k1) = b;                                                       % Store Parameter Estimates (Optional)
end
freq = 1./(2*diff(x0));                                                 % Calculate Frequencies By Cycle

figure(1)
plot(x, y)
hold on
plot(x0, zeros(size(x0)), '+r')
hold off
grid
axis([0  1    ylim])

The ‘freq’ calculation multiplies the difference by 2 because it is actually calculating the half-cycle frequencies otherwise. To calculate the full-cycle frequencies (the correct result), it is necessary to multiply the half-cycle frequencies by 2.

The plot simply demonstrates that the linear interpolation is approximately accurate. It is not necessary for the code.

I don’t have your data, but it should work with it. It might give an anomalous result for the end value because it may calculate a zero-crossing there even if there is not one. Just discard it. My code is reasonably robust but cannot account for absolutely all possible conditions with such real-world data.

Answer 2

[Greg Dionne]

Assuming you have the Signal Processing Toolbox and your data x is uniformly sampled with sample rate Fs:

z = hilbert(x);
f = angle(z(2:end).*conj(z(1:end-1))).*(Fs / (2*pi));

That should get you started.

If you're comfortable with filtering concepts, look up Hilbert filtering, you can often design a pair of filters that can filter out unwanted noise / signals all in one go.

Take a look at FMDEMOD if you are simply oscillating within a narrow frequency range.