diff options
Diffstat (limited to 'macros/intfilt.sci')
| -rw-r--r-- | macros/intfilt.sci | 196 |
1 files changed, 107 insertions, 89 deletions
diff --git a/macros/intfilt.sci b/macros/intfilt.sci index 15da95c..44fffbf 100644 --- a/macros/intfilt.sci +++ b/macros/intfilt.sci @@ -1,36 +1,54 @@ function [h, a]= intfilt(R, L, freqmult) - - + + // This function estimate Interpolated FIR Filter Design. // Calling Sequence // h=intfilt(R,L,freqmult) // [h a]=intfilt(R,L,freqmult) - + // Parameters // R: Samples. It should be numeric // L: bandlimited interpolation samples. It must be nonzero. - // freqmult: bandlimitedness of ALPHA times the Nyquist frequency, IT can be numeric or character ('B' or 'L', B is length + // freqmult: bandlimitedness of ALPHA times the Nyquist frequency, IT can be numeric or character ('B' or 'L', B is length // (N+1)*L-1 for N odd and (N+1)*L for N even) - - // h: linear phase FIR filter. - + + // h: linear phase FIR filter. + // Examples - // h=intfilt(20,10,'l') - // h=intfilt(20,10,1) - // + // h=intfilt(20,10,'l') // The output of this example has 220 columns ,so it is difficult to write it here. + // h=intfilt(20,10,1) // The output of this example has 220 columns ,so it is difficult to write it here. + + //h1=intfilt(2,3,'l'); + //OUTPUT : + // - 0.0625 0. 0.5625 1. 0.5625 0. - 0.0625 + + //h2=intfilt(4,1,1); + //OUTPUT : + // 0.3001054 0.6366198 0.9003163 1. 0.9003163 0.6366198 0.3001054 + // See also // Authors // Jitendra Singh - +funcprot(0); +[lhs,rhs]=argn(0); + +if (rhs~=3) then + error ("Wrong number of input arguments.") +end + +if (lhs<1 | lhs>2) then + error ("Wrong number of input arguments.") +end + if or(type(R)==10) | or(type(L)==10) then error ('Argument R and L must be numeric.') - + else - - - - + + + + if argn(2)==3 then if type(freqmult)==10 then typ=freqmult; @@ -39,51 +57,51 @@ function [h, a]= intfilt(R, L, freqmult) freqmult=double(freqmult); typ='b'; end - + end - + if freqmult==0 then h=repmat(%nan,[1,(2*R*L-1)]) a=1; else - - + + //typ(1)=='b' | typ(1)=='B' - + if convstr(typ(1), 'u') =='B' then n=2*R*L-1; - - + + if freqmult==1 then M=[R R 0 0]; F= [0 1/(2*R) 1/(2*R) 0.5]; else M=R*[1 1]; - - if type(freqmult)==10 then + + if type(freqmult)==10 then F=[0 98/2/R]; else F=[0 freqmult/2/R] end - + for f=(1/R):(1/R):.5, - + if type(freqmult)==10 then - F=[F f-(98/2/R) f+(98/2/R)]; + F=[F f-(98/2/R) f+(98/2/R)]; else - F=[F f-(freqmult/2/R) f+(freqmult/2/R)]; + F=[F f-(freqmult/2/R) f+(freqmult/2/R)]; end - + M=[M 0 0]; end; if (F(length(F))>.5), F(length(F))=.5; end; - end + end N=n-1; F=F*2; M=M - - + + if (max(F)>1) | (min(F)<0) error('Frequencies in F must be in range [0,1]') end @@ -107,9 +125,9 @@ end ftype = 0; differ = 0; -N = N+1; - -F=F(:)/2; M=M(:); W=sqrt(W(:)); +N = N+1; + +F=F(:)/2; M=M(:); W=sqrt(W(:)); dF = diff(F); if (length(F) ~= length(W)*2) @@ -118,9 +136,9 @@ end if or(dF<0), - + error('F frequency must be increasing') - + end @@ -142,67 +160,67 @@ Nodd = N-fix(N./2).*2; if ~Nodd - m=(0:L)+.5; + m=(0:L)+.5; else - m=(0:L); + m=(0:L); end - - + + k=m'; need_matrix = (~band) | (~weights); - - - - + + + + if need_matrix - - I1=k(:,ones(size(m,1),size(m,2)))+m(ones(size(k,1),size(k,2)),:); - I2=k(:,ones(size(m,1),size(m,2)))-m(ones(size(k,1),size(k,2)),:); + + I1=k(:,ones(size(m,1),size(m,2)))+m(ones(size(k,1),size(k,2)),:); + I2=k(:,ones(size(m,1),size(m,2)))-m(ones(size(k,1),size(k,2)),:); G=zeros(size(I1,1),size(I1,2)); end if Nodd k=k(2:length(k)); - b0=0; + b0=0; end; b=zeros(size(k,1),size(k,2)); - + dd=diff(F); - + if or(dd==0) & R==1 then - + h=repmat(%nan,[1,n]) a=1 - + else for s=1:2:length(F), - - m=(M(s+1)-M(s))/(F(s+1)-F(s)); - b1=M(s)-m*F(s); + + m=(M(s+1)-M(s))/(F(s+1)-F(s)); + b1=M(s)-m*F(s); if Nodd b0 = b0 + (b1*(F(s+1)-F(s)) + m/2*(F(s+1)*F(s+1)-F(s)*F(s)))* abs(W((s+1)/2)^2) ; end - + b=b(:) b = b+(m/(4*%pi*%pi)*(cos(2*%pi*k*F(s+1))-cos(2*%pi*k*F(s)))./(k.*k))* abs(W((s+1)/2)^2); - - - + + + b = b' + (F(s+1)*(m*F(s+1)+b1)*sinf(2*k*F(s+1))- F(s)*(m*F(s)+b1)*sinf(2*k*F(s)))* abs(W((s+1)/2)^2); if need_matrix - - + + mat=matrix((.5*F(s+1)*(sinf(2*I1*F(s+1))+sinf(2*I2*F(s+1)))- .5*F(s)*(sinf(2*I1*F(s))+sinf(2*I2*F(s))) ) * abs(W((s+1)/2)^2),size(G,1),size(G,2)) ; mat=mat'; - G=G+mat; - + G=G+mat; + end end; - + if Nodd b=[b0; b']; end; @@ -219,23 +237,23 @@ Nodd = N-fix(N./2).*2; h=[a(L+1:-1:2)/2; a(1); a(2:L+1)/2].'; else h=.5*[flipud(a); a].'; - end; - end; - - //typ(1)=='l' | typ(1)=='L' - + end; + end; + + //typ(1)=='l' | typ(1)=='L' + + - elseif convstr(typ(1), 'u') =='L' then - + if n==0 then h=ones(1,R) - return - end - + return + end + t=0:n*R+1; l=ones(n+1,length(t)); - + for i=1:n+1 for j=1:n+1 if (j~=i) then @@ -243,9 +261,9 @@ Nodd = N-fix(N./2).*2; end end end - + h=zeros(1,(n+1)*R); - + for i=0:R-1 for j=0:n h(j*R+i+1)=l((n-j)+1,round((n-1)/2*R+i+1)); @@ -260,15 +278,15 @@ if h(1) == 0, else error ('This type of filter is not recognized.') - - + + end a=1; - - + + end end - + endfunction @@ -276,16 +294,16 @@ endfunction ////// Supplementary function function y=sinf(x) - + for i=1:length(x) if x(i)==0 then y(i)=1; else - + y(i)=sin(%pi*x(i))/(%pi*x(i)); end - - end - + + end + y=y'; endfunction |