This file is indexed.

/usr/share/octave/packages/signal-1.2.2/specgram.m is in octave-signal 1.2.2-1build1.

This file is owned by root:root, with mode 0o644.

The actual contents of the file can be viewed below.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
## Copyright (C) 1999-2001 Paul Kienzle <pkienzle@users.sf.net>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

## usage: [S [, f [, t]]] = specgram(x [, n [, Fs [, window [, overlap]]]])
##
## Generate a spectrogram for the signal. This chops the signal into
## overlapping slices, windows each slice and applies a Fourier
## transform to determine the frequency components at that slice.
##
## x: vector of samples
## n: size of fourier transform window, or [] for default=256
## Fs: sample rate, or [] for default=2 Hz
## window: shape of the fourier transform window, or [] for default=hanning(n)
##    Note: window length can be specified instead, in which case
##    window=hanning(length)
## overlap: overlap with previous window, or [] for default=length(window)/2
##
## Return values
##    S is complex output of the FFT, one row per slice
##    f is the frequency indices corresponding to the rows of S.
##    t is the time indices corresponding to the columns of S.
##    If no return value is requested, the spectrogram is displayed instead.
##
## Example
##    x = chirp([0:0.001:2],0,2,500);  # freq. sweep from 0-500 over 2 sec.
##    Fs=1000;                  # sampled every 0.001 sec so rate is 1 kHz
##    step=ceil(20*Fs/1000);    # one spectral slice every 20 ms
##    window=ceil(100*Fs/1000); # 100 ms data window
##    specgram(x, 2^nextpow2(window), Fs, window, window-step);
##
##    ## 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] = specgram(x, fftn, Fs, window, window-step);
##    S = abs(S(2:fftn*4000/Fs,:)); # magnitude in range 0<f<=4000 Hz.
##    S = S/max(S(:));           # normalize magnitude so that max is 0 dB.
##    S = max(S, 10^(-40/10));   # clip below -40 dB.
##    S = min(S, 10^(-3/10));    # clip above -3 dB.
##    imagesc(t, f, flipud(log(S)));   # display in log scale
##
## 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,:))));

function [S_r, f_r, t_r] = specgram(x, n = min(256, length(x)), Fs = 2, window = hanning(n), overlap = ceil(length(window)/2))

  if nargin < 1 || nargin > 5
    print_usage;
  ## make sure x is a vector
  elseif columns(x) != 1 && rows(x) != 1
    error ("specgram data must be a vector");
  end
  if columns(x) != 1, x = x'; end

  ## if only the window length is given, generate hanning window
  if length(window) == 1, window = hanning(window); end

  ## should be extended to accept a vector of frequencies at which to
  ## evaluate the fourier transform (via filterbank or chirp
  ## z-transform)
  if length(n)>1, 
    error("specgram doesn't handle frequency vectors yet"); 
  endif

  ## compute window offsets
  win_size = length(window);
  if (win_size > n)
    n = win_size;
    warning ("specgram fft size adjusted to %d", n);
  end
  step = win_size - overlap;

  ## build matrix of windowed data slices
  offset = [ 1 : step : length(x)-win_size ];
  S = zeros (n, length(offset));
  for i=1:length(offset)
    S(1:win_size, i) = x(offset(i):offset(i)+win_size-1) .* window;
  endfor

  ## compute fourier transform
  S = fft (S);

  ## extract the positive frequency components
  if rem(n,2)==1
    ret_n = (n+1)/2;
  else
    ret_n = n/2;
  end
  S = S(1:ret_n, :);

  f = [0:ret_n-1]*Fs/n;
  t = offset/Fs;
  if nargout==0
    imagesc(t, f, 20*log10(abs(S)));
    set (gca (), "ydir", "normal");
    xlabel ("Time")
    ylabel ("Frequency")
  endif
  if nargout>0, S_r = S; endif
  if nargout>1, f_r = f; endif
  if nargout>2, t_r = t; endif

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] = specgram(x);

%! ## test of returned shape
%!assert (rows(S), 128)
%!assert (columns(f), rows(S))
%!assert (columns(t), columns(S))
%!test [S, f, t] = specgram(x');
%!assert (rows(S), 128)
%!assert (columns(f), rows(S));
%!assert (columns(t), columns(S));
%!error (isempty(specgram([])));
%!error (isempty(specgram([1, 2 ; 3, 4])));
%!error (specgram)

%!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] = specgram(x);
%! specgram(x, 2^nextpow2(window), Fs, window, window-step);
%! disp("shows a diagonal from bottom left to top right");
%! input("press enter:","s");
%!
%! ## test of returned values
%! S = specgram(x, 2^nextpow2(window), Fs, window, window-step);
%! imagesc(20*log10(flipud(abs(S))));
%! disp("same again, but this time using returned value");

%!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] = specgram(x, fftn, Fs, window, window-step);
%! S = abs(S(2:fftn*4000/Fs,:)); # magnitude in range 0<f<=4000 Hz.
%! S = S/max(max(S));         # normalize magnitude so that max is 0 dB.
%! S = max(S, 10^(-40/10));   # clip below -40 dB.
%! S = min(S, 10^(-3/10));    # clip above -3 dB.
%! imagesc(flipud(20*log10(S)));
%!
%! % The image contains a spectrogram of 'sample.wav'