This file is indexed.

/usr/share/octave/packages/signal-1.3.2/fir2.m is in octave-signal 1.3.2-5.

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
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
## Copyright (C) 2000 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/>.

## -*- texinfo -*-
## @deftypefn  {Function File} {@var{b} =} fir2 (@var{n}, @var{f}, @var{m})
## @deftypefnx {Function File} {@var{b} =} fir2 (@var{n}, @var{f}, @var{m}, @var{grid_n})
## @deftypefnx {Function File} {@var{b} =} fir2 (@var{n}, @var{f}, @var{m}, @var{grid_n}, @var{ramp_n})
## @deftypefnx {Function File} {@var{b} =} fir2 (@var{n}, @var{f}, @var{m}, @var{grid_n}, @var{ramp_n}, @var{window})
##
## Produce an order @var{n} FIR filter with arbitrary frequency response
## @var{m} over frequency bands @var{f}, returning the @var{n}+1 filter
## coefficients in @var{b}.  The vector @var{f} specifies the frequency band
## edges of the filter response and @var{m} specifies the magnitude response
## at each frequency.
##
## The vector @var{f} must be nondecreasing over the range [0,1], and the
## first and last elements must be 0 and 1, respectively.  A discontinuous
## jump in the frequency response can be specified by duplicating a band edge
## in @var{f} with different values in @var{m}.
##
## The resolution over which the frequency response is evaluated can be
## controlled with the @var{grid_n} argument.  The default is 512 or the
## next larger power of 2 greater than the filter length.
##
## The band transition width for discontinuities can be controlled with the
## @var{ramp_n} argument.  The default is @var{grid_n}/25.  Larger values
## will result in wider band transitions but better stopband rejection.
##
## An optional shaping @var{window} can be given as a vector with length
## @var{n}+1.  If not specified, a Hamming window of length @var{n}+1 is used.
##
## To apply the filter, use the return vector @var{b} with the @code{filter}
## function, for example @code{y = filter (b, 1, x)}.
##
## Example:
## @example
## f = [0, 0.3, 0.3, 0.6, 0.6, 1]; m = [0, 0, 1, 1/2, 0, 0];
## [h, w] = freqz (fir2 (100, f, m));
## plot (f, m, ";target response;", w/pi, abs (h), ";filter response;");
## @end example
## @seealso{filter, fir1}
## @end deftypefn

function b = fir2(n, f, m, grid_n, ramp_n, window)

  if nargin < 3 || nargin > 6
    print_usage;
  endif

  ## verify frequency and magnitude vectors are reasonable
  t = length(f);
  if t<2 || f(1)!=0 || f(t)!=1 || any(diff(f)<0)
    error ("fir2: frequency must be nondecreasing starting from 0 and ending at 1");
  elseif t != length(m)
    error ("fir2: frequency and magnitude vectors must be the same length");
  ## find the grid spacing and ramp width
  elseif (nargin>4 && length(grid_n)>1) || ...
         (nargin>5 && (length(grid_n)>1 || length(ramp_n)>1))
    error ("fir2: grid_n and ramp_n must be integers");
  endif
  if nargin < 4, grid_n=[]; endif
  if nargin < 5, ramp_n=[]; endif

  ## find the window parameter, or default to hamming
  w=[];
  if length(grid_n)>1, w=grid_n; grid_n=[]; endif
  if length(ramp_n)>1, w=ramp_n; ramp_n=[]; endif
  if nargin < 6, window=w; endif
  if isempty(window), window=hamming(n+1); endif
  if !isreal(window) || ischar(window), window=feval(window, n+1); endif
  if length(window) != n+1, error ("fir2: window must be of length n+1"); endif

  ## Default grid size is 512... unless n+1 >= 1024
  if isempty (grid_n)
    if n+1 < 1024
      grid_n = 512;
    else
      grid_n = n+1;
    endif
  endif

  ## ML behavior appears to always round the grid size up to a power of 2
  grid_n = 2 ^ nextpow2 (grid_n);

  ## Error out if the grid size is not big enough for the window
  if 2*grid_n < n+1
    error ("fir2: grid size must be greater than half the filter order");
  endif

  if isempty (ramp_n), ramp_n = fix (grid_n / 25); endif

  ## Apply ramps to discontinuities
  if (ramp_n > 0)
    ## remember original frequency points prior to applying ramps
    basef = f(:); basem = m(:);

    ## separate identical frequencies, but keep the midpoint
    idx = find (diff(f) == 0);
    f(idx) = f(idx) - ramp_n/grid_n/2;
    f(idx+1) = f(idx+1) + ramp_n/grid_n/2;
    f = [f(:);basef(idx)]';

    ## make sure the grid points stay monotonic in [0,1]
    f(f<0) = 0;
    f(f>1) = 1;
    f = unique([f(:);basef(idx)(:)]');

    ## preserve window shape even though f may have changed
    m = interp1(basef, basem, f);

    ## axis([-.1 1.1 -.1 1.1])
    ## plot(f,m,'-xb;ramped;',basef,basem,'-or;original;'); pause;
  endif

  ## interpolate between grid points
  grid = interp1(f,m,linspace(0,1,grid_n+1)');
  ## hold on; plot(linspace(0,1,grid_n+1),grid,'-+g;grid;'); hold off; pause;

  ## Transform frequency response into time response and
  ## center the response about n/2, truncating the excess
  if (rem(n,2) == 0)
    b = ifft([grid ; grid(grid_n:-1:2)]);
    mid = (n+1)/2;
    b = real ([ b([end-floor(mid)+1:end]) ; b(1:ceil(mid)) ]);
  else
    ## Add zeros to interpolate by 2, then pick the odd values below.
    b = ifft([grid ; zeros(grid_n*2,1) ;grid(grid_n:-1:2)]);
    b = 2 * real([ b([end-n+1:2:end]) ; b(2:2:(n+1))]);
  endif

  ## Multiplication in the time domain is convolution in frequency,
  ## so multiply by our window now to smooth the frequency response.
  ## Also, for matlab compatibility, we return return values in 1 row
  b = b(:)' .* window(:)';

endfunction

%% Test that the grid size is rounded up to the next power of 2
%!xtest
%! f = [0 0.6 0.6 1]; m = [1 1 0 0];
%! b9  = fir2 (30, f, m, 9);
%! b16 = fir2 (30, f, m, 16);
%! b17 = fir2 (30, f, m, 17);
%! b32 = fir2 (30, f, m, 32);
%! assert ( isequal (b9,  b16))
%! assert ( isequal (b17, b32))
%! assert (~isequal (b16, b17))

%% Test expected magnitudes of passbands, stopbands, and cutoff frequencies
%!test
%! f = [0, 0.7, 0.7, 1]; m = [0, 0, 1, 1];
%! b = fir2 (50, f, m);
%! h = abs (freqz (b, 1, [0, 0.7, 1], 2));
%! assert (h(1) <= 3e-3)
%! assert (h(2) <= 1/sqrt (2))
%! assert (h(3), 1, 2e-3)

%!test
%! f = [0, 0.25, 0.25, 0.75, 0.75, 1]; m = [0, 0, 1, 1, 0, 0];
%! b = fir2 (50, f, m);
%! h = abs (freqz (b, 1, [0, 0.25, 0.5, 0.75, 1], 2));
%! assert (h(1) <= 3e-3)
%! assert (h(2) <= 1/sqrt (2))
%! assert (h(3), 1, 2e-3)
%! assert (h(4) <= 1/sqrt (2))
%! assert (h(5) <= 3e-3)

%!test
%! f = [0, 0.45, 0.45, 0.55, 0.55, 1]; m = [1, 1, 0, 0, 1, 1];
%! b = fir2 (50, f, m);
%! h = abs (freqz (b, 1, [0, 0.45, 0.5, 0.55, 1], 2));
%! assert (h(1), 1, 2e-3)
%! assert (h(2) <= 1/sqrt (2))
%! assert (h(3) <= 1e-1)
%! assert (h(4) <= 1/sqrt (2))
%! assert (h(5), 1, 2e-3)

%!demo
%! f=[0, 0.3, 0.3, 0.6, 0.6, 1]; m=[0, 0, 1, 1/2, 0, 0];
%! [h, w] = freqz(fir2(100,f,m));
%! subplot(121);
%! plot(f,m,';target response;',w/pi,abs(h),';filter response;');
%! subplot(122);
%! plot(f,20*log10(m+1e-5),';target response (dB);',...
%!      w/pi,20*log10(abs(h)),';filter response (dB);');

%!demo
%! f=[0, 0.3, 0.3, 0.6, 0.6, 1]; m=[0, 0, 1, 1/2, 0, 0];
%! plot(f,20*log10(m+1e-5),';target response;');
%! hold on;
%! [h, w] = freqz(fir2(50,f,m,512,0));
%! plot(w/pi,20*log10(abs(h)),';filter response (ramp=0);');
%! [h, w] = freqz(fir2(50,f,m,512,25.6));
%! plot(w/pi,20*log10(abs(h)),';filter response (ramp=pi/20 rad);');
%! [h, w] = freqz(fir2(50,f,m,512,51.2));
%! plot(w/pi,20*log10(abs(h)),';filter response (ramp=pi/10 rad);');
%! hold off;

%!demo
%! % Classical Jakes spectrum
%! % X represents the normalized frequency from 0
%! % to the maximum Doppler frequency
%! asymptote = 2/3;
%! X = linspace(0,asymptote-0.0001,200);
%! Y = (1 - (X./asymptote).^2).^(-1/4);
%!
%! % The target frequency response is 0 after the asymptote
%! X = [X, asymptote, 1];
%! Y = [Y, 0, 0];
%!
%! plot(X,Y,'b;Target spectrum;');
%! hold on;
%! [H,F]=freqz(fir2(20, X, Y));
%! plot(F/pi,abs(H),'c;Synthesized spectrum (n=20);');
%! [H,F]=freqz(fir2(50, X, Y));
%! plot(F/pi,abs(H),'r;Synthesized spectrum (n=50);');
%! [H,F]=freqz(fir2(200, X, Y));
%! plot(F/pi,abs(H),'g;Synthesized spectrum (n=200);');
%! hold off;
%! title('Theoretical/Synthesized CLASS spectrum');
%! xlabel('Normalized frequency (Fs=2)');
%! ylabel('Magnitude');