summaryrefslogtreecommitdiff
path: root/macros/invfreqz.sci
diff options
context:
space:
mode:
Diffstat (limited to 'macros/invfreqz.sci')
-rw-r--r--macros/invfreqz.sci193
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));
+*/
generated by cgit (git 2.34.1) at 2025-03-23 23:28:06 +0000