Search code examples
matlabsignal-processingspectrogram

Matlab cos signal from defined frequency course

how come the spectrogram of this code is at its max at approximately 4400Hz instead of 2400Hz in the last timestep of the stft? (frequencyCourse(end) = 100, meshingOrder = 24 -> f = 2400)

startTime = 0; %s
endTime = 30; %s
startIAS = 15; %Hz
endIAS = 100; %Hz
meshingOrder = 24;
fs = 100000; %Hz
t = startTime:1/fs:endTime-1/fs;
frequencyCourse = linspace(startIAS, endIAS, length(t));
signal = cos(2*pi*meshingOrder*frequencyCourse.*t);
spectrogram(signal, hanning(2^13), 0, 2^14, fs, 'yaxis')

Here's a picture:

enter image description here

It works fine as long as I use chirp instead of my self constructed signal, it's not an option though, since more specific courses are to come.

Solution

  • Summary

    The problem is that the instantaneous phase is the integral of the instantaneaous frequency with respect to time, not the instantaneous frequency multiplied by time.

    You should compute the signal as

    signal = cos(2*pi*meshingOrder*cumtrapz(t, frequencyCourse));

    What your code does

    In your example, you seem to want to generate a linear chirp with initial frequency meshingOrder*startIAS and final frequency meshingOrder*endIAS. But that's not the code is doing.

    In your computed signal, the instantaneous phase is the argument to the cos function:

    2*pi*meshingOrder*frequencyCourse.*t

    Since the variable frequencyCourse increases from meshingOrder*startIAS at start time (which is 0) to meshingOrder*endIAS at end time, this can be expressed as

    2*pi*(A+B*t).*t

    where A = meshingOrder*startIAS and B = meshingOrder*(endIAS-startIAS)/endTime. Differentiating the instantanous phase with respect to t gives an instantaneous frequency

    A + 2*B*t

    that is

    meshingOrder*startIAS + 2*meshingOrder*(endIAS-startIAS)/endTime * t

    As you can see, the problem is the factor 2 here. At end time the instantaneous frequency is

    meshingOrder*startIAS + 2*meshingOrder*(endIAS-startIAS)

    that is

    2*meshingOrder*endIAS - meshingOrder*startIAS

    In your example this is 4440 Hz, which is in accordance with your observed value.

    What the code should do

    For a linear chirp (or a chirp with any other simple frequency variation, such as a quadratic or exponential) you could compute the correct instantaneous phase that gives rise to the desired instantaneous frequency. See for example here. This is also what the chirp function internally does.

    But you seem to be want to deal with arbitrary frequency courses. To do that, given arbitrary t, just compute the argument of the cos as the cumulative integral of frequencyCourse with respect to t. This is easily done with cumtrapz:

    signal = cos(2*pi*meshingOrder*cumtrapz(t, frequencyCourse));

    Changing this line in your example gives the following figure, which has the expected frequency variation from 360 Hz to 2400 Hz:

    enter image description here