diff options
Diffstat (limited to 'macros/invfreqs.sci')
| -rw-r--r-- | macros/invfreqs.sci | 202 |
1 files changed, 108 insertions, 94 deletions
diff --git a/macros/invfreqs.sci b/macros/invfreqs.sci index 32e732b..99c721b 100644 --- a/macros/invfreqs.sci +++ b/macros/invfreqs.sci @@ -1,96 +1,110 @@ -function [B,A,C] = invfreqs(H,F,nB,nA,W,iter,tol,trace) -//Fit filter B(s)/A(s)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] = invfreqs(H,F,nB,nA,W,iter,tol,trace) -//[B,A,C] = invfreqs(H,F,nB,nA,W) -//[B,A,C] = invfreqs(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(s)/A(s)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 = [1/2 1]; -//A = [1 1]; -//w = linspace(0,4,128); -//H = freqs(B,A,w); -//[Bh,Ah, C] = invfreqs(H,w,1,1); -//Bh = -// -// 0.50000 1.00000 -// -//Ah = -// -// 1.0000 1.0000 -// -//C = -3.0964e-15 - -funcprot(0); -lhs = argn(1) -rhs = argn(2) -if (rhs < 4 | rhs > 8) -error("Wrong number of input arguments.") -end - -select(rhs) - - case 4 then - if(lhs==1) - B = callOctave("invfreqs",H,F,nB,nA) - elseif(lhs==2) - [B, A] = callOctave("invfreqs",H,F,nB,nA) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqs",H,F,nB,nA) - else - error("Wrong number of output argments.") - end - - case 5 then - if(lhs==1) - B = callOctave("invfreqs",H,F,nB,nA,W) - elseif(lhs==2) - [B, A] = callOctave("invfreqs",H,F,nB,nA,W) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqs",H,F,nB,nA,W) - else - error("Wrong number of output argments.") - end - case 6 then - if(lhs==1) - B = callOctave("invfreqs",H,F,nB,nA,W,iter) - elseif(lhs==2) - [B, A] = callOctave("invfreqs",H,F,nB,nA,W,iter) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqs",H,F,nB,nA,W,iter) - else - error("Wrong number of output argments.") - end - case 7 then - if(lhs==1) - B = callOctave("invfreqs",H,F,nB,nA,W,iter,tol) - elseif(lhs==2) - [B, A] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol) - else - error("Wrong number of output argments.") - end - case 8 then - if(lhs==1) - B = callOctave("invfreqs",H,F,nB,nA,W,iter,tol,trace) - elseif(lhs==2) - [B, A] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol,trace) - elseif(lhs==3) - [B, A, C] = callOctave("invfreqs",H,F,nB,nA,W,iter,tol,trace) - else - error("Wrong number of output argments.") - end - end +// 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] = invfreqs(H,F,nB,nA) + :[B,A] = invfreqs(H,F,nB,nA,W) + :[B,A] = invfreqs(H,F,nB,nA,W,iter,tol,'trace') + + Fit filter B(s)/A(s)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 + + H: desired complex frequency response + + F: frequency (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] = invfreqs(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, 's', varargin); + endfunction +/* +demo + B = [1 0 0]; + A = [1 6 15 15]/15; + w = linspace(0, 8, 128); + [H0 ,_ ] = freqz(B, A, w); + Nn = (rand(size(w,1),size(w,2),'normal')+%i*rand(size(w,1),size(w,2),'normal'))/sqrt(2); + order = length(A) - 1; + [Bh, Ah, Sig0] = invfreqs(H0, w, [length(B)-1 2], length(A)-1); + [Hh ,_ ] = freqz(Bh,Ah,w); + [BLS, ALS, SigLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "LS"); + [HLS,_] = freqz(BLS, ALS, w); + [BTLS, ATLS, SigTLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "TLS"); + [HTLS,_] = freqz(BTLS, ATLS, w); + [BMLS, AMLS, SigMLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "QR"); + [HMLS,_] = freqz(BMLS, AMLS, w); + plot(w,[abs(H0); abs(Hh)]) + xlabel("Frequency (rad/sec)"); + ylabel("Magnitude"); + legend('Original','Measured'); + err = norm(H0-Hh); + disp(sprintf('L2 norm of frequency response error = %f',err)); +*/ +/* Octave version + B = [1 0 0]; + A = [1 6 15 15]/15; + w = linspace(0, 8, 128); + [H0 ,_ ] = freqz(B, A, w); + Nn = (randn(size(w,1),size(w,2))+i*randn(size(w,1),size(w,2)))/sqrt(2); + order = length(A) - 1; + [Bh, Ah, Sig0] = invfreqs(H0, w, [length(B)-1 2], length(A)-1); + [Hh ,_ ] = freqz(Bh,Ah,w); + [BLS, ALS, SigLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "LS"); + [HLS,_] = freqz(BLS, ALS, w); + [BTLS, ATLS, SigTLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "TLS"); + [HTLS,_] = freqz(BTLS, ATLS, w); + [BMLS, AMLS, SigMLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [], "method", "QR"); + [HMLS,_] = freqz(BMLS, AMLS, w); + plot(w,[abs(H0); abs(Hh)]) + xlabel("Frequency (rad/sec)"); + ylabel("Magnitude"); + legend('Original','Measured'); + err = norm(H0-Hh); + disp(sprintf('L2 norm of frequency response error = %f',err)); +*/ + + + |