Clear Filters
Clear Filters

Why are the results of my spectrogram sinusoidal?

4 views (last 30 days)
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;
Any help would be extremely gratefully received!

Accepted Answer

J. Webster
J. Webster on 21 Apr 2016
Edited: J. Webster on 21 Apr 2016
You don't want just the real part of the fft, you want the magnitude of it.
s_real = real(s)
s_mag = abs(s)
and use that in your plot.
J. Webster
J. Webster on 21 Apr 2016
Edited: J. Webster on 21 Apr 2016
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

Sign in to comment.

More Answers (1)

Star Strider
Star Strider on 21 Apr 2016
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:
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)
freqz(sos, 4028, fs); % Filter Bode Plot
yf = filtfilt(sos,g, y); % Filter Signal (Phase-Neutral)
env = abs(hilbert(yf)); % Calculate Envelope
plot(t, y)
title('Unfiltered Signal')
axis([0 0.2 ylim]) % Axis Resolution (Comment-Out To See Entire Signal)
plot(t, yf)
hold on
plot(t, env, '-r', t, -env, '-r')
hold off
title('Bandpass-Filtered Signal')
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.
Star Strider
Star Strider on 21 Apr 2016
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 on 21 Apr 2016
I will try the filter bank method. Strange that the spectrogram appears steady in figure 2, but could be the windowing.

Sign in to comment.


Community Treasure Hunt

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

Start Hunting!