Spectrogram behaves unexpectedly when changing input values systematically

12 ビュー (過去 30 日間)
Mike Denny
Mike Denny 2017 年 6 月 26 日
コメント済み: Mike Denny 2017 年 7 月 3 日
I hope this will be fun for you to investigate and play with. Weird stuff happens. If weird stuff doesn't happen for you, then maybe my version of MATLAB (R2013a) is to blame.
I want to make a spectrogram of a noisy signal. Due to the signal's complexity, I've been messing around with a matlab example to understand how the inputs affect the output and what I expect to happen based on the help info vs. what I see actually plotted. This is the example (I added the 'yaxis' bit):
t=0:0.001:2; % 2 secs @ 1kHz sample rate
y=chirp(t,100,1,200,'q'); % Start @ 100Hz, cross 200Hz at t=1sec
spectrogram(y,128,120,128,1E3,'yaxis'); % Display the spectrogram
title('Quadratic Chirp: start at 100Hz and cross 200Hz at t=1sec');
I can tell, then, that they're using the form "spectrogram(X,WINDOW,NOVERLAP,NFFT,Fs,FREQLOCATION)."
So my understanding is that the 2nd argument is the number of time-axis segments to split the data into, 3rd is not super clear actually but the default is 50%. 50% of, I initially assumed, the adjacent two segments. The later statement of "NOVERLAP is the number of samples each segment of X overlaps." Ok, well in the default case of 8 segments and 50% of the samples, how can you make 8 segments when each segment is supposed to have 50% of the samples in it?? And why does it show 7 segments? This will show you what I mean:
spectrogram(y,[],[],128,1E3,'yaxis');
Then lets decide we want 128 segments like they did, and we'll just be ok with the default NOVERLAP. Well, plot that and you'll see 29 segments, which isn't 50% of anything I can see.
Now put NOVERLAP to 20, and we get 17 segments. Change it to 5 and you'll get 15! Plot the original again, then lets look at what happens when we double the number of segments and bump the NOVERLAP a bit:
spectrogram(y,256,200,128,1E3,'yaxis');
That doesn't look right! Maybe what I thought about NFFT was unclear. Lets double that...nope...at least that seems to behave properly. That's the number of segments the sample frequency spectrum is broken in to. Just for fun you can defy logic one last time with this:
spectrogram(y,256,200,50,1E3,'yaxis');
How am I supposed to select the appropriate values for this function? Is there something wrong with it? For my real case I'll have a sampling frequency of 72 kHz and will have applied a zero-phase Butterworth low pass filter with cutoff at about 7 kHz.

採用された回答

Shane L
Shane L 2017 年 6 月 26 日
The 2nd argument (if it is a single integer), WINDOW, is not the number of time segments into which to split the data, but the length of each time segment. To see this, try, for example, a modified version of the example code:
t = 0:0.001:2-0.001; % 2 secs @ 1kHz sample rate
y = chirp(t,100,1,200,'q'); % start @ 100Hz, cross 200Hz at t = 1 sec
spectrogram(y,200,0,200,1e3,'yaxis'); % display the spectrogram
The signal has a length of 2000 samples. We set the second argument, WINDOW, to 200, which specifies that each segment contains 200 samples. We also set the 3rd input, NOVERLAP, to zero for this example. Since the signal with 2000 samples is split into segments of 200 samples, there are 10 time segments in the resulting spectrogram, as shown in the figure below.
The 3rd argument, NOVERLAP, specifies the number of samples that overlap between two adjacent segments in the calculation of the short-time Fourier transform. In the example above, we set NOVERLAP to zero, so no samples were overlapped. By default, 50% of the samples are overlapped. To see this, try, for example:
t = 0:0.001:2-0.001; % 2 secs @ 1kHz sample rate
y = chirp(t,100,1,200,'q'); % start @ 100Hz, cross 200Hz at t = 1 sec
spectrogram(y,200,[],200,1e3,'yaxis'); % display the spectrogram
Now, you will see 19 time segments in the resulting spectrogram, as shown in the figure below. The reason is because of the overlap. The first segment includes samples 1-200; the second segment includes samples 101-300 (50% overlap of 200 samples is 100 samples); the third segment includes samples 201-400; up to the 19th segment, which includes samples 1801-2000.
I don't know if this was part of your question, but the 4th argument, NFFT, is the number of points used in the calculation of the FFT for each segment. If NFFT is less than WINDOW, then not all of the points in each segment are used in the FFT. If NFFT is greater than WINDOW, then the segment is zero-padded before computing the FFT. Try changing the value of NFFT (for example, 10 vs. 1000) and plotting the spectrogram, and you will see the frequency resolution change, as shown in the two figures below.
For more information, you can consult the MATLAB help page for spectrogram ; several examples demonstrate how to compute the number of resulting time and frequency bins based on WINDOW, NOVERLAP, and NFFT.
  4 件のコメント
Shane L
Shane L 2017 年 6 月 27 日
Hi Mike, I got a chance to access R2014b. MATLAB introduced additional code into the R2015a version of spectrogram that shifts the time series by a half-interval and adds an additional point at the end, which explains the differences in the plots. See below for the block of code in R2015a:
% Surf encodes data points at the edges and takes the color from
% the last edge so we need to add an additional point so that surf
% does the right thing. This is important especially when
% spectrogram has only one estimate (e.g. window length = signal
% length). Although this issue also exists on the frequency
% direction we will not add an extra row since we will never
% encounter a single frequency point estimate. For the plot we set
% time values to be at: nwin/2-a/2, nwin/2+a/2, nwin/2+3a/2,
% nwin/2+5a/2 ... where a is the number of new samples for each
% segment (i.e. nwind-noverlap). For the case of zero overlap this
% corresponds to 0, nwind, 2*nwind, ...
a = nwind - noverlap;
t = [(nwind/2-a/2)/Fs t+((a/2)/Fs)];
Pxx = [Pxx Pxx(:,end)];
If you would like your plots to be in style of R2015a and later, then edit spectrogram and insert that block of code in between lines 195 and 196 (if the code in your version is the same as in R2014b, otherwise, find the code as shown below:
[Pxx,W] = compute_PSD(win,y,nfft,f,options,Fs,esttype); %#ok<NASGU>
%%%INSERT LINES HERE
displayspectrogram(t,f,double(Pxx),isFsnormalized,faxisloc);
I tested my first example code using R2014b and it now matches R2015a. I hope this helps!
Mike Denny
Mike Denny 2017 年 7 月 3 日
Hi Shane,
I didn't get a notification that you also commented so sorry for the delay. But yep, the Matlab revision was definitely the issue. I didn't try the code, however, since I had the ability to upgrade to 2015b.
Thanks for the help.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeTime-Frequency Analysis についてさらに検索

製品

Community Treasure Hunt

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

Start Hunting!

Translated by