Why are the results of my spectrogram sinusoidal?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Harvey
am 21 Apr. 2016
Kommentiert: Harvey
am 22 Apr. 2016
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/162099/image.png)
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):
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/162102/image.png)
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!
0 Kommentare
Akzeptierte Antwort
J. Webster
am 21 Apr. 2016
Bearbeitet: J. Webster
am 21 Apr. 2016
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.
3 Kommentare
J. Webster
am 21 Apr. 2016
Bearbeitet: J. Webster
am 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
Weitere Antworten (1)
Star Strider
am 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:
[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.
3 Kommentare
Star Strider
am 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.
Siehe auch
Kategorien
Mehr zu Spectral Measurements finden Sie in Help Center und File Exchange
Produkte
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!