Why are the results of my spectrogram sinusoidal?

I am trying to determine the amplitude envelope of specific frequencies over time, from a sample of an instrument (a trumpet). I use the spectrogram function to find the amplitude of each frequency bin over time. Knowing that a harmonic lies within a certain bin, I can plot that particular bin against time to produce the following figure:

Although I can make out the general envelope, as marked in red, this method is much less accurate than what I was expecting.

I can't understand why the results are sinusoidal. Is this a result of aliasing, or something else entirely?

If we compare this plot, to a plot of all frequencies, we notice that there is no sinusoidal fluctuation in the power of the particular harmonic (The harmonic plotted in the first figure has been boxed):

So where is difference in the values used to plot the consistent power we see here, and the values that produce the sinusoidal plot?

For reference I have attached the audio file I am using. Bellow are the functions I used to get these results. Row 13 refers to the bin containing the specific harmonic;

 [y,fs]=audioread('Trumpet_C4.wav');
 [s,f,t]=spectrogram(y,hann(2048),1024,2048,fs,'yaxis');
 s_real=real(s);
 plot(t,s_real(13,:))
 figure(2)
 spectrogram(y,hann(2048),1024,2048,fs,'taxi's')

Any help would be extremely gratefully received!

Answer 1

[J. Webster]

You don't want just the real part of the fft, you want the magnitude of it.

change

s_real = real(s)

to

s_mag = abs(s)

and use that in your plot.

[Harvey]

I think that did it! Thanks.

[J. Webster]

I believe that what's happening is that when you call spectrogram with no output arguments, it plots the PSD of the signal, but when you call it with output arguments, it returns the complex sffts of the signal.

Note that the abs(s) is not equal to the PSD, it just has the same behavior. If you want the PSD instead of just the magnitude, I think you can just use

PSD = 1/(fs*2048) .* s_mag.^2    %2048 is the number of points in the fft

[Harvey]

Thank you!

Answer 2

[Star Strider]

I am not certain what you want to do. Just filtering the fundamental frequency gives it as a function of time, and the envelope is relatively easy to calculate:

[y,fs]=audioread('Trumpet_C4.wav');
t = linspace(0, length(y), length(y))/fs;                       % Time Vector

Fn = fs/2;                                                      % Nyquist Frequency

Wp = [0.2  0.4]*1E+3/Fn;                                        % Filter Passband
Ws = [0.5  2.0].*Wp;                                            % Filter Stopband
Rp = 10;                                                        % Passband Ripple
Rs = 45;                                                        % Stopband Ripple
[n,Wn] = buttord(Wp, Ws, Rp, Rs);                               % Filter Order
[b,a] = butter(n,Wn);                                           % Transfer Function Coefficient Vectors
[sos,g] = tf2sos(b,a);                                          % Second-Order-Section (For Stability)

figure(1)
freqz(sos, 4028, fs);                                           % Filter Bode Plot

yf = filtfilt(sos,g, y);                                        % Filter Signal (Phase-Neutral)

env = abs(hilbert(yf));                                         % Calculate Envelope

figure(2)
subplot(2,1,1)
plot(t, y)
title('Unfiltered Signal')
grid
axis([0  0.2    ylim])                                          % Axis Resolution (Comment-Out To See Entire Signal)
subplot(2,1,2)
plot(t, yf)
hold on
plot(t, env, '-r',  t, -env, '-r')
hold off
title('Bandpass-Filtered Signal')
grid
axis([0  0.2    ylim])                                          % Axis Resolution (Comment-Out To See Entire Signal)
xlabel('Time (s)')

It’s a signal with a fundamental frequency and (apparently 13) harmonics, so you would have to filter each one, and calculate the envelope for each one. They seem to be evenly-spaced, so you could determine the spacing and then design a bank of filters (see How to shape the spectrum of an audio waveform? to filter each frequency.

[Harvey]

Thank you very much. I will try this method, although I can't see it being easily applied to multiple harmonics. My final goal, is a Power against Time plot of the first 64 harmonics (there are only 13 in that figure as it has been cropped). I will be applying this to multiple instrument samples.

Do you have any suggestions as to why the method I used returned the results it did?

Why does the spectrogram function plot a steady Power over time (figure 2) if the first figure shows sinusoidal output or vice versa?

The ideal outcome is to plot the Power values used in figure 2 against their time values.

Thanks

[Star Strider]

My pleasure.

There are several ways to get the individual harmonics. I would use the filter-bank method to get them all as individual records. I’ve had success with it in other applications, and it probably produces the most accurate results for what you want to do.

I went through the spectrogram documentation and couldn’t figure out the reason it gave you the results it did. I suspect it has to do with windowing, overlap and other necessary operations it has to do to calculate its results.

[Harvey]

I will try the filter bank method. Strange that the spectrogram appears steady in figure 2, but could be the windowing.

Thanks