octave_packages/signal-1.1.3/fir2.m

   1 ## Copyright (C) 2000 Paul Kienzle <pkienzle@users.sf.net>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## usage: b = fir2(n, f, m [, grid_n [, ramp_n]] [, window])
  17 ##
  18 ## Produce an FIR filter of order n with arbitrary frequency response,
  19 ## returning the n+1 filter coefficients in b.
  20 ##
  21 ## n: order of the filter (1 less than the length of the filter)
  22 ## f: frequency at band edges
  23 ##    f is a vector of nondecreasing elements in [0,1]
  24 ##    the first element must be 0 and the last element must be 1
  25 ##    if elements are identical, it indicates a jump in freq. response
  26 ## m: magnitude at band edges
  27 ##    m is a vector of length(f)
  28 ## grid_n: length of ideal frequency response function
  29 ##    defaults to 512, should be a power of 2 bigger than n
  30 ## ramp_n: transition width for jumps in filter response
  31 ##    defaults to grid_n/20; a wider ramp gives wider transitions
  32 ##    but has better stopband characteristics.
  33 ## window: smoothing window
  34 ##    defaults to hamming(n+1) row vector
  35 ##    returned filter is the same shape as the smoothing window
  36 ##
  37 ## To apply the filter, use the return vector b:
  38 ##       y=filter(b,1,x);
  39 ## Note that plot(f,m) shows target response.
  40 ##
  41 ## Example:
  42 ##   f=[0, 0.3, 0.3, 0.6, 0.6, 1]; m=[0, 0, 1, 1/2, 0, 0];
  43 ##   [h, w] = freqz(fir2(100,f,m));
  44 ##   plot(f,m,';target response;',w/pi,abs(h),';filter response;');
  45
  46 function b = fir2(n, f, m, grid_n, ramp_n, window)
  47
  48   if nargin < 3 || nargin > 6
  49     print_usage;
  50   endif
  51
  52   ## verify frequency and magnitude vectors are reasonable
  53   t = length(f);
  54   if t<2 || f(1)!=0 || f(t)!=1 || any(diff(f)<0)
  55     usage("frequency must be nondecreasing starting from 0 and ending at 1");
  56   endif
  57   if t != length(m)
  58     usage("frequency and magnitude vectors must be the same length");
  59   endif
  60
  61   ## find the grid spacing and ramp width
  62   if (nargin>4 && length(grid_n)>1) || \
  63         (nargin>5 && (length(grid_n)>1 || length(ramp_n)>1))
  64     usage("grid_n and ramp_n must be integers");
  65   endif
  66   if nargin < 4, grid_n=512; endif
  67   if nargin < 5, ramp_n=grid_n/20; endif
  68
  69   ## find the window parameter, or default to hamming
  70   w=[];
  71   if length(grid_n)>1, w=grid_n; grid_n=512; endif
  72   if length(ramp_n)>1, w=ramp_n; ramp_n=grid_n/20; endif
  73   if nargin < 6, window=w; endif
  74   if isempty(window), window=hamming(n+1); endif
  75   if !isreal(window) || ischar(window), window=feval(window, n+1); endif
  76   if length(window) != n+1, usage("window must be of length n+1"); endif
  77
  78   ## make sure grid is big enough for the window
  79   if 2*grid_n < n+1, grid_n = 2^nextpow2(n+1); endif
  80
  81   ## Apply ramps to discontinuities
  82   if (ramp_n > 0)
  83     ## remember original frequency points prior to applying ramps
  84     basef = f(:); basem = m(:);
  85
  86     ## separate identical frequencies, but keep the midpoint
  87     idx = find (diff(f) == 0);
  88     f(idx) = f(idx) - ramp_n/grid_n/2;
  89     f(idx+1) = f(idx+1) + ramp_n/grid_n/2;
  90     f = [f(:);basef(idx)]';
  91
  92     ## make sure the grid points stay monotonic in [0,1]
  93     f(f<0) = 0;
  94     f(f>1) = 1;
  95     f = unique([f(:);basef(idx)(:)]');
  96
  97     ## preserve window shape even though f may have changed
  98     m = interp1(basef, basem, f);
  99
 100     # axis([-.1 1.1 -.1 1.1])
 101     # plot(f,m,'-xb;ramped;',basef,basem,'-or;original;'); pause;
 102   endif
 103
 104   ## interpolate between grid points
 105   grid = interp1(f,m,linspace(0,1,grid_n+1)');
 106   # hold on; plot(linspace(0,1,grid_n+1),grid,'-+g;grid;'); hold off; pause;
 107
 108   ## Transform frequency response into time response and
 109   ## center the response about n/2, truncating the excess
 110   if (rem(n,2) == 0)
 111     b = ifft([grid ; grid(grid_n:-1:2)]);
 112     mid = (n+1)/2;
 113     b = real ([ b([end-floor(mid)+1:end]) ; b(1:ceil(mid)) ]);
 114   else
 115     ## Add zeros to interpolate by 2, then pick the odd values below.
 116     b = ifft([grid ; zeros(grid_n*2,1) ;grid(grid_n:-1:2)]);
 117     b = 2 * real([ b([end-n+1:2:end]) ; b(2:2:(n+1))]);
 118   endif
 119   ## Multiplication in the time domain is convolution in frequency,
 120   ## so multiply by our window now to smooth the frequency response.
 121   if rows(window) > 1
 122     b = b .* window;
 123   else
 124     b = b' .* window;
 125   endif
 126 endfunction
 127
 128 %!demo
 129 %! f=[0, 0.3, 0.3, 0.6, 0.6, 1]; m=[0, 0, 1, 1/2, 0, 0];
 130 %! [h, w] = freqz(fir2(100,f,m));
 131 %! subplot(121);
 132 %! plot(f,m,';target response;',w/pi,abs(h),';filter response;');
 133 %! subplot(122);
 134 %! plot(f,20*log10(m+1e-5),';target response (dB);',...
 135 %!      w/pi,20*log10(abs(h)),';filter response (dB);');
 136
 137 %!demo
 138 %! f=[0, 0.3, 0.3, 0.6, 0.6, 1]; m=[0, 0, 1, 1/2, 0, 0];
 139 %! plot(f,20*log10(m+1e-5),';target response;');
 140 %! hold on;
 141 %! [h, w] = freqz(fir2(50,f,m,512,0));
 142 %! plot(w/pi,20*log10(abs(h)),';filter response (ramp=0);');
 143 %! [h, w] = freqz(fir2(50,f,m,512,25.6));
 144 %! plot(w/pi,20*log10(abs(h)),';filter response (ramp=pi/20 rad);');
 145 %! [h, w] = freqz(fir2(50,f,m,512,51.2));
 146 %! plot(w/pi,20*log10(abs(h)),';filter response (ramp=pi/10 rad);');
 147 %! hold off;
 148
 149 %!demo
 150 %! % Classical Jakes spectrum
 151 %! % X represents the normalized frequency from 0
 152 %! % to the maximum Doppler frequency
 153 %! asymptote = 2/3;
 154 %! X = linspace(0,asymptote-0.0001,200);
 155 %! Y = (1 - (X./asymptote).^2).^(-1/4);
 156 %!
 157 %! % The target frequency response is 0 after the asymptote
 158 %! X = [X, asymptote, 1];
 159 %! Y = [Y, 0, 0];
 160 %!
 161 %! title('Theoretical/Synthesized CLASS spectrum');
 162 %! xlabel('Normalized frequency (Fs=2)');
 163 %! ylabel('Magnitude');
 164 %!
 165 %! plot(X,Y,'b;Target spectrum;');
 166 %! hold on;
 167 %! [H,F]=freqz(fir2(20, X, Y));
 168 %! plot(F/pi,abs(H),'c;Synthesized spectrum (n=20);');
 169 %! [H,F]=freqz(fir2(50, X, Y));
 170 %! plot(F/pi,abs(H),'r;Synthesized spectrum (n=50);');
 171 %! [H,F]=freqz(fir2(200, X, Y));
 172 %! plot(F/pi,abs(H),'g;Synthesized spectrum (n=200);');
 173 %! hold off;
 174 %! xlabel(''); ylabel(''); title('');