## -*- texinfo -*- ## @deftypefn {} {@var{s} =} spectrogram (@var{x}) ## @deftypefnx {} {@var{s} =} spectrogram (@var{x}, @var{window}) ## @deftypefnx {} {@var{s} =} spectrogram (@var{x}, @var{window}, @var{noverlap}) ## @deftypefnx {} {@var{s} =} spectrogram (@var{x}, @var{window}, @var{noverlap}, @var{nfft}) ## @deftypefnx {} {[@var{s}, @var{w}, @var{t}] =} spectrogram (@dots{}) ## @deftypefnx {} {[@var{s}, @var{f}, @var{t}] =} spectrogram (@dots{}, @var{fs}) ## @deftypefnx {} {[@var{s}, @var{w}, @var{t}] =} spectrogram (@var{x}, @var{window}, @var{noverlap}, @var{w}) ## @deftypefnx {} {[@var{s}, @var{f}, @var{t}] =} spectrogram (@var{x}, @var{window}, @var{noverlap}, @var{f}, @var{fs}) ## @deftypefnx {} {[@dots{}, @var{ps}] =} spectrogram (@dots{}) ## @deftypefnx {} {[@dots{}] =} spectrogram (@dots{}, "reassigned"}) ## @deftypefnx {} {[@dots{}, @var{ps}, @var{fc}, @var{tc}] =} spectrogram (@dots{}) ## @deftypefnx {} {[@dots{}] =} spectrogram (@dots{}, @var{freqrange}) ## @deftypefnx {} {[@dots{}] =} spectrogram (@dots{}, @var{spectrumtype}) ## @deftypefnx {} {[@dots{}] =} spectrogram (@dots{}, "minthreshold", @var{threshold}) ## @deftypefnx {} {} spectrogram (@dots{}) ## @deftypefnx {} {} spectrogram (@dots{}, @var{freqloc}) ## Compute a spectrogram for the signal @var{x}. ## ## All arguments except @var{x} are optional. The data is divided into ## segments. If @var{window} is a vector, each segment has the same ## length as @var{window} and is multiplied by @var{window} before ## (optional) zero-padding and calculation of its periodogram. If ## @var{window} is a scalar, each segment has a length of @var{window} ## and a Hamming window is used. ## ## If @var{noverlap} is an integer scalar, it is the number of ## elements of overlap for each segment. If it is missing or empty, ## then an overlap of 50% is used for each segment. ## ## If both @var{window} and @var{noverlap} are missing or empty, then ## @var{x} is divided into 8 segments that overlap by 50%. ## ## The spectral density is the mean of the periodograms, scaled so that ## area under the spectrum is the same as the mean square of the ## data. This equivalence is supposed to be exact, but in practice there ## is a mismatch of up to 0.5% when comparing area under a periodogram ## with the mean square of the data. ## ## The optional argument @var{nfft} must be a non-negative integer ## value specifying the number of points to use when computing the FFT ## of each segment. By default, @var{nfft} is computed to be the ## @code{max (256, 2 ^ nextpow2 (seg_len))}. ## ## If the fifth argument is a positive scalar it is taken to be the ## sampling frequency. Otherwise a value of @math{2\pi} is used. ## ## FIXME: Need explanation of optional input (and output) arguments ## @var{w} and @var{f}. ## ## FIXME: Need explanation of output arguments @var{ps}, @var{fc} and ## @var{tc}. ## ## The optional argument @var{freqrange} may be one of the following: ## ## @table @qcode ## @item "onesided" ## The frequency range of the spectrum is zero up to but not including ## @code{fs/2}. Power from negative frequencies is added to the ## positive side of the spectrum, but not at zero or Nyquist (@code{fs/2}) ## frequencies. This keeps power equal in time and spectral domains. ## See reference [2]. ## ## @item "twosided" ## The frequency range of the spectrum is @code{-fs/2} to @code{fs/2}, ## with negative frequencies stored in "wrap around" order after the positive ## frequencies; e.g. frequencies for a 10-point 'twosided' ## spectrum are 0 0.1 0.2 0.3 0.4 0.5 -0.4 -0.3 -0.2 -0.1 ## ## @item "centered" ## The same as "twosided" but with the first half of the spectrum ## swapped with second half to put the zero-frequency value in the ## middle. @seealso{fftshift}. ## @end table ## ## If @var{x} is real, the default for @var{freqrange} is ## @qcode{"onesided"}, otherwise the default @qcode{"twosided"}. ## ## The optional argument @var{spectrumtype} may be one of the ## following ## ## @table @qcode ## @item "psd" ## @item "power" ## @end @table ## ## The optional argument @var{freqloc} may be either @qcode{"xaxis"}, ## to display frequency on the x-axis (the default) when plotting ## results, or @qcode{"yaxis"} to display it on the y-axis. ## ## If the optional name/value pair @code{"minthreshold"} and a ## corresponding scalar value, FIXME: DO WHAT? ## ## If no return values are requested, the results are plotted and ## no values are returned. ## ## FIXME: Does the following info belong here? If so, does it also ## belong in the pwelch doc string? ## ## The choice of window defines the time-frequency resolution. In ## speech for example, a wide window shows more harmonic detail while a ## narrow window averages over the harmonic detail and shows more ## formant structure. The shape of the window is not so critical so long ## as it goes gradually to zero on the ends. ## ## Step size (which is window length minus overlap) controls the ## horizontal scale of the spectrogram. Decrease it to stretch, or ## increase it to compress. Increasing step size will reduce time ## resolution, but decreasing it will not improve it much beyond the ## limits imposed by the window size (you do gain a little bit, ## depending on the shape of your window, as the peak of the window ## slides over peaks in the signal energy). The range 1-5 msec is good ## for speech. ## ## FFT length controls the vertical scale. Selecting an FFT length ## greater than the window length does not add any information to the ## spectrum, but it is a good way to interpolate between frequency ## points which can make for prettier spectrograms. ## ## After you have generated the spectral slices, there are a number of ## decisions for displaying them. First the phase information is ## discarded and the energy normalized: ## ## S = abs(S); S = S/max(S(:)); ## ## Then the dynamic range of the signal is chosen. Since information in ## speech is well above the noise floor, it makes sense to eliminate any ## dynamic range at the bottom end. This is done by taking the max of ## the magnitude and some minimum energy such as minE=-40dB. Similarly, ## there is not much information in the very top of the range, so ## clipping to a maximum energy such as maxE=-3dB makes sense: ## ## S = max(S, 10^(minE/10)); S = min(S, 10^(maxE/10)); ## ## The frequency range of the FFT is from 0 to the Nyquist frequency of ## one half the sampling rate. If the signal of interest is band ## limited, you do not need to display the entire frequency range. In ## speech for example, most of the signal is below 4 kHz, so there is no ## reason to display up to the Nyquist frequency of 10 kHz for a 20 kHz ## sampling rate. In this case you will want to keep only the first 40% ## of the rows of the returned S and f. More generally, to display the ## frequency range [minF, maxF], you could use the following row index: ## ## idx = (f >= minF & f <= maxF); ## ## Then there is the choice of colormap. A brightness varying colormap ## such as copper or bone gives good shape to the ridges and valleys. A ## hue varying colormap such as jet or hsv gives an indication of the ## steepness of the slopes. The final spectrogram is displayed in log ## energy scale and by convention has low frequencies on the bottom of ## the image: ## ## imagesc(t, f, flipud(log(S(idx,:)))); ## @end deftypefn function [s, varargout] = spectrogram (x, varargin) if (nargin < 1) print_usage (); endif ## Force x to be column vector. if (rows (x) == 1) x = x(:); endif [window, noverlap, nfft, fs, f, w, info] = parse_args (x, varargin{:}); x_len = info.x_len; seg_len = info.seg_len; compute_f = false; compute_w = false; if (nargout > 2 || nargout == 0) if (info.have_f_arg || info.have_fs_arg) compute_f = true; else compute_w = true; endif endif compute_ps = (nargout == 0 || nargout > 3); compute_fc = (nargout > 4); compute_tc = (nargout > 4); f_val = []; t_val = []; w_val = w; ps_val = []; fc_val = []; tc_val = []; ## Calculate results. step = seg_len - noverlap; offset = 1:step:x_len-seg_len+1; n_ffts = numel (offset); tmp_s = zeros (nfft, n_ffts); if (strcmp (info.freqrange, "onesided")) if (! rem (nfft, 2)) # one-sided, nfft is even psd_len = nfft / 2 + 1; else # one-sided, nfft is odd psd_len = (nfft + 1) / 2; endif else # two-sided (and shifted) psd_len = nfft; endif ## If pwelch were fixed so that it worked column-wise on a matrix ## argument, then we could store the windowed values of X in a ## single large matrix and then call fft and pwelch each just one ## time. But this efficiencly would also require (possibly much) ## more memory. s_val = zeros (psd_len, n_ffts); if (compute_ps) ps_val = zeros (psd_len, n_ffts); endif if (compute_fc) fc_val = zeros (psd_len, n_ffts); endif if (compute_tc) tc_val = zeros (psd_len, n_ffts); endif for i = 1:n_ffts start_seg = offset(i); end_seg = offset(i)+seg_len-1; xx = x(start_seg:end_seg); tmp_s = fft (xx .* window, nfft); s_val(:,i) = tmp_s(1:psd_len,:); if (compute_ps) ps_val(:,i) = pwelch (xx, window, 0, nfft, info.freqrange, info.spectrumtype); endif ## FIXME: compute FC_VAL, and TC_VAL if needed. Is this the right ## place to do it? How? if (compute_fc) fc_val(:,i) = 1; endif if (compute_tc) tc_val(:,i) = 1; endif endfor if (compute_ps && info.have_minthreshold_arg) ps_val (10*log10 (ps_val) <= info.threshold) = 0; endif ## Frequency. freq = (fs / nfft) * (0:psd_len-1)'; ## FIXME: is this correct? What is W supposed to be? if (compute_w) if (info.have_w_arg) w_val = w; else w_val = freq; endif endif if (compute_f) if (info.have_f_arg) f_val = f; else f_val = freq; endif endif ## Centers of the segments: t_val = (seg_len/2 + (0:n_ffts-1)*(seg_len-noverlap)); if (info.have_fs_arg) t_val = t_val / fs; endif if (nargout > 0) s = s_val; if (compute_f) varargout{1} = f_val; varargout{2} = t_val; elseif (compute_w) varargout{1} = w_val; varargout{2} = t_val; endif if (compute_ps) varargout{3} = ps_val; endif if (compute_fc) varargout{4} = fc_val; endif if (compute_tc) varargout{5} = tc_val; endif else ## Although the primary computation is the short time Fourier ## transform, the default plot is the PSD for each segment, shown ## in decibels. db_val = 10 * log10 (ps_val+eps); if (info.have_fs_arg) ## We can display units of time since we know the sample rate. sample_axis_label = "Time (sec)"; freq_axis_label = "Frequency (Hz)"; sample_lim = [t_val(1), t_val(end)]; freq_lim = f_val; else sample_axis_label = "Samples"; freq_axis_label = "Normalized Frequency"; sample_lim = [t_val(1), t_val(end)]; freq_lim = [0, 1]; endif if (strcmp (info.freqloc, "xaxis")) db_val = db_val'; axdir = "xdir"; xlab = freq_axis_label; ylab = sample_axis_label; xlim = freq_lim; ylim = sample_lim; else axdir = "ydir"; xlab = sample_axis_label; ylab = freq_axis_label; xlim = sample_lim; ylim = freq_lim; endif imagesc (xlim, ylim, db_val); set (gca (), "ydir", "normal", "xlabel", xlab, "ylabel", ylab); colorbar ("ylabel", "Power/Frequency (db/Hz)"); endif endfunction ## The rules for determining argument meanings and default parameter ## and option values for SPECTROGRAM are as complicated as just about any ## function I've ever encountered. function [window, noverlap, nfft, fs, f, w, info] = parse_args (x, varargin) x_len = length (x); ## Find the first option argument in varargin. nvargs = numel (varargin); non_option_args = nargin; for i = 1:nvargs if (ischar (varargin{i})); ## Including X. non_option_args = i; break; endif endfor if (non_option_args > 5) evalin ("caller", "print_usage"); endif info.have_fs_arg = false; info.have_f_arg = false; info.have_w_arg = false; ## FS must be the fifth argument. if (non_option_args == 5) fs = varargin{4}; if (! (isscalar (fs) && fs > 0)) error ("spectrogram: FS must be a scalar greater than 0"); endif info.have_fs_arg = true; else fs = 2 * pi; endif if (non_option_args > 1) window = varargin{1}; else window = []; endif if (! (isempty (window) || (isscalar (window) && window > 0 && round (window) == window) || isvector (window))) error ("spectrogram: WINDOW must be a vector or or an integer greater than 0"); endif if (non_option_args > 2) noverlap = varargin{2}; else noverlap = []; endif if (! (isempty (noverlap) || (isscalar (noverlap) && noverlap >= 0 && round (noverlap) == noverlap))) error ("spectrogram: NOVERLAP must be an integer greater than or equal to 0"); endif ## At this point, WINDOW is empty, a positive scalar, or a vector, ## and NOVERLAP is empty, or a non-negative scalar. if (isempty (window)) if (isempty (noverlap)) ## 8 segments, 50% overlap. window = fix (x_len * 2 / 9); else ## Choose window such that there are 8 segments that overlap by ## noverlap elements. window = fix ((x_len + 7 * noverlap) / 8); endif endif if (isscalar (window)) seg_len = window; elseif (isvector (window)) seg_len = numel (window); else error ("spectrogram: internal error"); endif if (isempty (noverlap)) noverlap = fix (seg_len / 2); endif if (noverlap >= seg_len) error ("spectrogram: NOVERLAP must be less than window length"); endif if (isscalar (window)) ## Make Hamming window. xx = seg_len - 1; window = 0.54 - 0.46 * cos ((2 * pi / xx) * (0:xx).'); endif nfft = []; f = NaN (2, 1); w = NaN (2, 1); if (non_option_args > 3) argval = varargin{3}; if (isempty (argval)) ## Use default value for NFFT. elseif (isscalar (argval)) ## Scalar value must be NFFT. nfft = argval; if (! (nfft > 0 && round (nfft) == nfft)) error ("spectrogram: NFFT must be an integer greater than 0"); endif elseif (isvector (argval)) if (info.have_fs_arg) f = argval; info.have_f_arg = true; else w = argval; info.have_w_arg = true; endif else error ("spectrogram: invalid value for NFFT, F, or W") endif endif if (non_option_args < 4 || isempty (nfft)) nfft = max (256, 2 ^ nextpow2 (seg_len)); endif info.have_freqrange_arg = false; info.have_spectrumtype_arg = false; info.have_trace_arg = false; info.have_minthreshold_arg = false; info.have_freqloc_arg = false; info.reassigned = false; if (iscomplex (x)) info.freqrange = "twosided"; else info.freqrange = "onesided"; endif info.spectrumtype = "psd"; info.trace = "mean"; info.threshold = -Inf; info.freqloc = "xaxis"; vargidx = non_option_args; while (vargidx <= nvargs) argval = varargin{vargidx}; switch (tolower (argval)) case "reassigned" info.reassigned = true; case {"onesided", "twosided", "centered"} if (info.have_freqrange_arg) error ("spectrogram: multiple freqrange options specified"); endif info.have_freqrange_arg = true; info.freqrange = argval; case {"psd", "power"} if (info.have_spectrumtype_arg) error ("spectrogram: multiple spectrumtype options specified"); endif info.have_spectrumtype_arg = true; info.spectrumtype = argval; case {"minhold", "maxhold"} if (info.have_trace_arg) error ("spectrogram: multiple trace options specified"); endif info.have_trace_arg = true; info.trace = argval; case "minthreshold" if (info.have_minthreshold_arg) error ("spectrogram: multiple minthreshold options specified"); endif info.have_minthreshold_arg = true; vargidx++; if (vargidx <= nvargs) info.threshold = varargin{vargidx}; if (! isscalar (info.threshold)) error ("spectrogram: minthreshold must be a scalar value"); endif else error ("spectrogram: expecting value to follow 'minthreshold' argument"); endif case {"xaxis", "yaxis"} if (info.have_freqloc_arg) error ("spectrogram: multiple freqloc options specified"); endif info.have_freqloc_arg = true; info.freqloc = argval; otherwise error ("spectrogram: unrecognized option argument"); endswitch vargidx++; endwhile info.x_len = x_len; info.seg_len = seg_len; endfunction %!shared S,f,t,x %! Fs=1000; %! x = chirp([0:1/Fs:2],0,2,500); # freq. sweep from 0-500 over 2 sec. %! step=ceil(20*Fs/1000); # one spectral slice every 20 ms %! window=ceil(100*Fs/1000); # 100 ms data window %! [S, f, t] = spectrogram(x); %! ## test of returned shape %!assert (rows(S), 128) %!assert (columns(f), rows(S)) %!assert (columns(t), columns(S)) %!test [S, f, t] = spectrogram(x'); %!assert (rows(S), 128) %!assert (columns(f), rows(S)); %!assert (columns(t), columns(S)); %!error (isempty(spectrogram([]))); %!error (isempty(spectrogram([1, 2 ; 3, 4]))); %!error (spectrogram) %!demo %! Fs=1000; %! x = chirp([0:1/Fs:2],0,2,500); # freq. sweep from 0-500 over 2 sec. %! step=ceil(20*Fs/1000); # one spectral slice every 20 ms %! window=ceil(100*Fs/1000); # 100 ms data window %! %! ## test of automatic plot %! [S, f, t] = spectrogram(x); %! spectrogram(x, window, window-step, 2^nextpow2(window), Fs); %!#demo # FIXME: Enable once we have an audio file to demo %! ## Speech spectrogram %! [x, Fs] = auload(file_in_loadpath("sample.wav")); # audio file %! step = fix(5*Fs/1000); # one spectral slice every 5 ms %! window = fix(40*Fs/1000); # 40 ms data window %! fftn = 2^nextpow2(window); # next highest power of 2 %! [S, f, t] = spectrogram(x, window, window-step, fftn, Fs); %! S = abs(S(2:fftn*4000/Fs,:)); # magnitude in range 0