diff options
Diffstat (limited to 'macros/invfreqz.sci')
| -rw-r--r-- | macros/invfreqz.sci | 193 |
1 files changed, 103 insertions, 90 deletions
diff --git a/macros/invfreqz.sci b/macros/invfreqz.sci index d051499..433100b 100644 --- a/macros/invfreqz.sci +++ b/macros/invfreqz.sci @@ -1,94 +1,107 @@ -function [B,A,C] = invfreqz(H,F,nB,nA,W,iter,tol,trace) -//Fit filter B(z)/A(z)to the complex frequency response H at frequency points F. A and B are real polynomial coefficients of order nA and nB. -//Calling Sequence -//[B,A,C] = invfreqz(H,F,nB,nA,W,iter,tol,trace) -//[B,A,C] = invfreqz(H,F,nB,nA,W) -//[B,A,C] = invfreqz(H,F,nB,nA) -//Parameters -//H: desired complex frequency response. -//F: frequency (must be same length as H). -//nB: order of the numerator polynomial B. -//nA: order of the denominator polynomial A. -//W: vector of weights (must be same length as F). -//Description -//This is an Octave function. -//Fit filter B(z)/A(z)to the complex frequency response H at frequency points F. A and B are real polynomial coefficients of order nA and nB. -//Optionally, the fit-errors can be weighted vs frequency according to the weights W. -//Note: all the guts are in invfreq.m -//Examples -//[B,A] = butter(4,1/4); -//[H,F] = freqz(B,A); -//[Bh,Ah,C] = invfreq(H,F,4,4) -//Bh = -// -// 0.010209 0.040838 0.061257 0.040838 0.010209 -// -//Ah = -// -// 1.00000 -1.96843 1.73586 -0.72447 0.12039 -// -//C = -7.7065e-15 +// Copyright (C) 2018 - IIT Bombay - FOSSEE +// This file must be used under the terms of the CeCILL. +// This source file is licensed as described in the file COPYING, which +// you should have received as part of this distribution. The terms +// are also available at +// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt +// Original Source : https://octave.sourceforge.io/ +// Modifieded by: Abinash Singh Under FOSSEE Internship +// Last Modified on : 3 Feb 2024 +// Organization: FOSSEE, IIT Bombay +// Email: toolbox@scilab.in +// FIXME: check invfreq.sci for todo's +/* + : [B,A] = invfreqz(H,F,nB,nA) ¶ + : [B,A] = invfreqz(H,F,nB,nA,W) ¶ + : [B,A] = invfreqz(H,F,nB,nA,W,iter,tol,'trace') ¶ -funcprot(0); -lhs = argn(1) -rhs = argn(2) -if (rhs < 4 | rhs > 8) -error("Wrong number of input arguments.") -end + Fit filter B(z)/A(z)to the complex frequency response H at frequency points F. -select(rhs) - - case 4 then - if(lhs==1) - B = callOctave("invfreqz",H,F,nB,nA) - elseif(lhs==2) - [B, A] = callOctave("invfreqz",H,F,nB,nA) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqz",H,F,nB,nA) - else - error("Wrong number of output argments.") - end + A and B are real polynomial coefficients of order nA and nB. Optionally, the fit-errors can be weighted vs frequency according to the weights W. - case 5 then - if(lhs==1) - B = callOctave("invfreqz",H,F,nB,nA,W) - elseif(lhs==2) - [B, A] = callOctave("invfreqz",H,F,nB,nA,W) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqz",H,F,nB,nA,W) - else - error("Wrong number of output argments.") - end - case 6 then - if(lhs==1) - B = callOctave("invfreqz",H,F,nB,nA,W,iter) - elseif(lhs==2) - [B, A] = callOctave("invfreqz",H,F,nB,nA,W,iter) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqz",H,F,nB,nA,W,iter) - else - error("Wrong number of output argments.") - end - case 7 then - if(lhs==1) - B = callOctave("invfreqz",H,F,nB,nA,W,iter,tol) - elseif(lhs==2) - [B, A] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol) - else - error("Wrong number of output argments.") - end - case 8 then - if(lhs==1) - B = callOctave("invfreqz",H,F,nB,nA,W,iter,tol,trace) - elseif(lhs==2) - [B, A] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol,trace) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqz",H,F,nB,nA,W,iter,tol,trace) - else - error("Wrong number of output argments.") - end - end -endfunction + Note: all the guts are in invfreq.m + + H: desired complex frequency response + + F: normalized frequency (0 to pi) (must be same length as H) + + nA: order of the denominator polynomial A + + nB: order of the numerator polynomial B + W: vector of weights (must be same length as F) + +*/ +// Dependencies +// invfreq +function [B, A, SigN] = invfreqz(H, F, nB, nA, W, iter, tol, tr, varargin) + + if nargin < 9 + varargin = {}; + if nargin < 8 + tr = ''; + if nargin < 7 + tol = []; + if nargin < 6 + iter = []; + if nargin < 5 + W = ones(1,length(F)); + end + end + end + end + end + // now for the real work + [B, A, SigN] = invfreq(H, F, nB, nA, W, iter, tol, tr, 'z', varargin); + +endfunction +/* +demo +order = 9; //order of test filter +//going to 10 or above leads to numerical instabilities and large errors +fc = 1/2; // sampling rate / 4 +n = 128; // frequency grid size +// butterworth filter of order 9 and fc=0.5 +B0 = [5.1819e-03 4.6637e-02 1.8655e-01 4.3528e-01 6.5292e-01 6.5292e-01 4.3528e-01 1.8655e-01 4.6637e-02 5.1819e-03]; +A0 = [ 1.0000e+00 -8.6736e-16 1.2010e+00 -7.7041e-16 4.0850e-01 -1.7013e-16 4.2661e-02 -9.0155e-18 9.6666e-04 -5.3661e-20]; +[H0, w] = freqz(B0, A0, n); +Nn = (rand(size(w,1),size(w,2),'normal')+%i*rand(size(w,1),size(w,2),'normal'))/sqrt(2); +[Bh, Ah, Sig0] = invfreqz(H0, w, order, order); +[Hh, wh] = freqz(Bh, Ah, n); +[BLS, ALS, SigLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "LS"); +[HLS _ ] = freqz(BLS, ALS, n); +[BTLS, ATLS, SigTLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "TLS"); +[HTLS _ ]= freqz(BTLS, ATLS, n); +[BMLS, AMLS, SigMLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "QR"); +[HMLS _ ] = freqz(BMLS, AMLS, n); +plot(w,[abs(H0) abs(Hh)]) +xlabel("Frequency (rad/sample)"); +ylabel("Magnitude"); +legend('Original','Measured'); +err = norm(H0-Hh); +disp(sprintf('L2 norm of frequency response error = %f',err)); + +*/ +/* +order = 9; +fc = 1/2; +n = 128; +B0 = [5.1819e-03 4.6637e-02 1.8655e-01 4.3528e-01 6.5292e-01 6.5292e-01 4.3528e-01 1.8655e-01 4.6637e-02 5.1819e-03]; +A0 = [ 1.0000e+00 -8.6736e-16 1.2010e+00 -7.7041e-16 4.0850e-01 -1.7013e-16 4.2661e-02 -9.0155e-18 9.6666e-04 -5.3661e-20]; +[H0, w] = freqz(B0, A0, n); +Nn = (randn(size(w,1),size(w,2))+i*randn(size(w,1),size(w,2)))/sqrt(2); +[Bh, Ah, Sig0] = invfreqz(H0, w, order, order); +[Hh, wh] = freqz(Bh, Ah, n); +[BLS, ALS, SigLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "LS"); +[HLS _ ] = freqz(BLS, ALS, n); +[BTLS, ATLS, SigTLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "TLS"); +[HTLS _ ]= freqz(BTLS, ATLS, n); +[BMLS, AMLS, SigMLS] = invfreqz(H0+1e-5*Nn, w, order, order, [], [], [], [], "method", "QR"); +[HMLS _ ] = freqz(BMLS, AMLS, n); +plot(w,[abs(H0) abs(Hh)]) +xlabel("Frequency (rad/sample)"); +ylabel("Magnitude"); +legend('Original','Measured'); +err = norm(H0-Hh); +disp(sprintf('L2 norm of frequency response error = %f',err)); +*/ |