1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
|
// 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));
*/
|
