In reply to Soren message, I enclose the patchset for the three files,
invfreq*.m. I tried to keep the interface constant. The new output value is
an approximation of the noise level, evaluated from the misfit between
simulated and observed H(s).
Regards
Pascal
--
Save money, use Scarlet Internet solutions: as from EUR 25/month... visit
http://www.scarlet.be/
--- invfreq.m 2010-09-28 10:10:29.000000000 +0200
+++ /users/dupuis/matlab/octave-forge/signal-1.0.11/invfreq.m 2010-10-06
10:27:19.000000000 +0200
@@ -61,27 +61,63 @@
% *final debugging
% 2007-08-03 Rolf Schirmacher <[email protected]>
% *replaced == by strcmp() for character string comparison
+% 2010-10-04 Pascal Dupuis <[email protected]>
+% *added the ability to specify zeroes in B
% TODO: implement Steiglitz-McBride iterations
% TODO: improve numerical stability for high order filters (matlab is a bit
better)
% TODO: modify to accept more argument configurations
-function [B,A] = invfreq(H,F,nB,nA,W,iter,tol,tr,plane)
- n = max(nA,nB);
+function [B, A, SigN] = invfreq(H, F, nB, nA, W, iter, tol, tr, plane,
varargin)
+ if length(nB) > 1, zB = nB(2); nB = nB(1); else zB = 0; end
+ n = max(nA, nB);
m = n+1; mA = nA+1; mB = nB+1;
nF = length(F);
if nF ~= length(H), disp('invfreqz: length of H and F must be the same');
end;
- if nargin < 5, W = ones(1,nF); end;
+ if nargin < 5 || isempty(W), W = ones(1, nF); end;
if nargin < 6, iter = []; end
if nargin < 7 tol = []; end
- if nargin < 8, tr = ''; end
+ if nargin < 8 || isempty(tr), tr = ''; end
if nargin < 9, plane = 'z'; end
+ if nargin < 10, varargin = {}; end
if iter~=[], disp('no implementation for iter yet'),end
if tol ~=[], disp('no implementation for tol yet'),end
if plane ~= 'z' & plane ~= 's', disp('invfreqz: Error in plane argument'),
end
- Ruu = zeros(mB,mB); Ryy = zeros(nA,nA); Ryu = zeros(nA,mB);
- Pu = zeros(mB,1); Py = zeros(nA,1);
+ [reg, prop ] = parseparams(varargin);
+ %# should we normalise freqs to avoid matrices with rank deficiency ?
+ norm = false;
+ %# by default, use Ordinary Least Square to solve normal equations
+ method = 'LS';
+ if length(prop) > 0
+ indi = 1; while indi <= length(prop)
+ switch prop{indi}
+ case 'norm'
+ if indi < length(prop) && ~ischar(prop{indi+1}),
+ norm = logical(prop{indi+1});
+ prop(indi:indi+1) = [];
+ continue
+ else
+ norm = true; prop(indi) = [];
+ continue
+ end
+ case 'method'
+ if indi < length(prop) && ischar(prop{indi+1}),
+ method = prop{indi+1};
+ prop(indi:indi+1) = [];
+ continue
+ else
+ error('invfreq.m: incorrect/missing method argument');
+ end
+ otherwise %# FIXME: just skip it for now
+ disp(sprintf("Ignoring unkown argument %s", varargin{indi}));
+ indi = indi + 1;
+ end
+ end
+ end
+
+ Ruu = zeros(mB, mB); Ryy = zeros(nA, nA); Ryu = zeros(nA, mB);
+ Pu = zeros(mB, 1); Py = zeros(nA,1);
if strcmp(tr,'trace')
disp(' ')
disp('Computing nonuniformly sampled, equation-error, rational filter.');
@@ -90,48 +126,118 @@
end
s = sqrt(-1)*F;
- if strcmp(plane,'z')
- if max(F)>pi || min(F)<0
+ switch plane
+ case 'z'
+ if max(F) > pi || min(F) < 0
disp('hey, you frequency is outside the range 0 to pi, making my own')
- F = linspace(0,pi,length(H));
+ F = linspace(0, pi, length(H));
s = sqrt(-1)*F;
end
s = exp(-s);
- end;
+ case 's'
+ if max(F) > 1e6 && n > 5,
+ if ~norm,
+ disp('Be carefull, there are risks of generating singular matrices');
+ disp('Call invfreqs as (..., "norm", true) to avoid it');
+ else
+ Fmax = max(F); s = sqrt(-1)*F/Fmax;
+ end
+ end
+ end
- for k=1:nF
+ for k=1:nF,
Zk = (s(k).^[0:n]).';
Hk = H(k);
aHks = Hk*conj(Hk);
Rk = (W(k)*Zk)*Zk';
rRk = real(Rk);
- Ruu = Ruu+rRk(1:mB,1:mB);
- Ryy = Ryy+aHks*rRk(2:mA,2:mA);
- Ryu = Ryu+real(Hk*Rk(2:mA,1:mB));
- Pu = Pu+W(k)*real(conj(Hk)*Zk(1:mB));
- Py = Py+(W(k)*aHks)*real(Zk(2:mA));
+ Ruu = Ruu + rRk(1:mB, 1:mB);
+ Ryy = Ryy + aHks*rRk(2:mA, 2:mA);
+ Ryu = Ryu + real(Hk*Rk(2:mA, 1:mB));
+ Pu = Pu + W(k)*real(conj(Hk)*Zk(1:mB));
+ Py = Py + (W(k)*aHks)*real(Zk(2:mA));
end;
+ Rr = ones(length(s), mB+nA); Zk = s;
+ for k = 1:min(nA, nB),
+ Rr(:, 1+k) = Zk;
+ Rr(:, mB+k) = -Zk.*H;
+ Zk = Zk.*s;
+ end
+ for k = 1+min(nA, nB):max(nA, nB)-1,
+ if k <= nB, Rr(:, 1+k) = Zk; end
+ if k <= nA, Rr(:, mB+k) = -Zk.*H; end
+ Zk = Zk.*s;
+ end
+ k = k+1;
+ if k <= nB, Rr(:, 1+k) = Zk; end
+ if k <= nA, Rr(:, mB+k) = -Zk.*H; end
+
+ %# complex to real equation system -- this ensures real solution
+ Rr = Rr(:, 1+zB:end);
+ Rr = [real(Rr); imag(Rr)]; Pr = [real(H(:)); imag(H(:))];
+ %# normal equations -- keep for ref
+ %# Rn= [Ruu(1+zB:mB, 1+zB:mB), -Ryu(:, 1+zB:mB)'; -Ryu(:, 1+zB:mB), Ryy];
+ %# Pn= [Pu(1+zB:mB); -Py];
+
+ switch method
+ case {'ls' 'LS'}
+ %# avoid scaling errors with Theta = R\P;
+ %# [Q, R] = qr([Rn Pn]); Theta = R(1:end, 1:end-1)\R(1:end, end);
+ [Q, R] = qr([Rr Pr], 0); Theta = R(1:end-1, 1:end-1)\R(1:end-1, end);
+ %# SigN = R(end, end-1);
+ SigN = R(end, end);
+ case {'tls' 'TLS'}
+ % [U, S, V] = svd([Rn Pn]);
+ % SigN = S(end, end-1);
+ % Theta = -V(1:end-1, end)/V(end, end);
+ [U, S, V] = svd([Rr Pr], 0);
+ SigN = S(end, end);
+ Theta = -V(1:end-1, end)/V(end, end);
+ case {'mls' 'MLS' 'qr' 'QR'}
+ % [Q, R] = qr([Rn Pn], 0);
+ %# solve the noised part -- DO NOT USE ECONOMY SIZE !
+ % [U, S, V] = svd(R(nA+1:end, nA+1:end));
+ % SigN = S(end, end-1);
+ % Theta = -V(1:end-1, end)/V(end, end);
+ %# unnoised part -- remove B contribution and back-substitute
+ % Theta = [R(1:nA, 1:nA)\(R(1:nA, end) - R(1:nA, nA+1:end-1)*Theta)
+ % Theta];
+ %# solve the noised part -- economy size OK as #rows > #columns
+ [Q, R] = qr([Rr Pr], 0);
+ eB = mB-zB; sA = eB+1;
+ [U, S, V] = svd(R(sA:end, sA:end));
+ %# noised (A) coefficients
+ Theta = -V(1:end-1, end)/V(end, end);
+ %# unnoised (B) part -- remove A contribution and back-substitute
+ Theta = [R(1:eB, 1:eB)\(R(1:eB, end) - R(1:eB, sA:end-1)*Theta)
+ Theta];
+ SigN = S(end, end);
+ otherwise
+ error("invfreq: unknown method %s", method);
+ end
- R = [Ruu,-Ryu';-Ryu,Ryy];
- P = [Pu;-Py];
- Theta = R\P;
+ B = [zeros(zB, 1); Theta(1:mB-zB)].';
+ A = [1; Theta(mB-zB+(1:nA))].';
- B = Theta(1:mB)';
- A = [1 Theta(mB+1:mB+nA)'];
if strcmp(plane,'s')
B = B(mB:-1:1);
A = A(mA:-1:1);
+ if norm, %# Frequencies were normalised -- unscale coefficients
+ Zk = Fmax.^[n:-1:0].';
+ for k = nB:-1:1+zB, B(k) = B(k)/Zk(k); end
+ for k = nA:-1:1, A(k) = A(k)/Zk(k); end
+ end
end
endfunction
%!demo
-%! order = 12; % order of test filter
+%! order = 6; % order of test filter
%! fc = 1/2; % sampling rate / 4
%! n = 128; % frequency grid size
-%! [B,A] = butter(order,fc);
-%! [H,w] = freqz(B,A,n);
-%! [Bh,Ah] = invfreq(H,w,order,order);
-%! [Hh,wh] = freqz(Bh,Ah,n);
+%! [B, A] = butter(order,fc);
+%! [H, w] = freqz(B,A,n);
+%! [Bh, Ah] = invfreq(H,w,order,order);
+%! [Hh, wh] = freqz(Bh,Ah,n);
%! xlabel("Frequency (rad/sample)");
%! ylabel("Magnitude");
%! plot(w,[abs(H), abs(Hh)])
--- invfreqs.m 2010-08-05 18:19:51.000000000 +0200
+++ /users/dupuis/matlab/octave-forge/signal-1.0.11/invfreqs.m 2010-10-06
09:33:23.000000000 +0200
@@ -49,9 +49,11 @@
% TODO: check invfreq.m for todo's
-function [B,A] = invfreqs(H,F,nB,nA,W,iter,tol,tr)
+function [B, A, SigN] = invfreqs(H,F,nB,nA,W,iter,tol,tr, varargin)
-if nargin < 8
+if nargin < 9
+ varargin = {};
+ if nargin < 8
tr = '';
if nargin < 7
tol = [];
@@ -62,22 +64,35 @@
end
end
end
+ end
end
% now for the real work
-[B,A] = invfreq(H,F,nB,nA,W,iter,tol,tr,'s');
+[B, A, SigN] = invfreq(H, F,nB, nA, W, iter, tol, tr, 's', varargin{:});
endfunction
%!demo
%! B = [1/2 1];
+%! B = [1 0 0];
%! A = [1 1];
-%! w = linspace(0,4,128);
-%! H = freqs(B,A,w);
-%! [Bh,Ah] = invfreqs(H,w,1,1);
+%! %#A = [1 36 630 6930 51975 270270 945945 2027025 2027025]/2027025;
+%! A = [1 21 210 1260 4725 10395 10395]/10395;
+%! A = [1 6 15 15]/15;
+%! w = linspace(0, 8, 128);
+%! H0 = freqs(B, A, w);
+%! Nn = (randn(size(w))+j*randn(size(w)))/sqrt(2);
+%! order = length(A) - 1;
+%! [Bh, Ah, Sig0] = invfreqs(H0, w, [length(B)-1 2], length(A)-1);
%! Hh = freqs(Bh,Ah,w);
+%! [BLS, ALS, SigLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [], [],
"method", "LS");
+%! HLS = freqs(BLS, ALS, w);
+%! [BTLS, ATLS, SigTLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [],
[], "method", "TLS");
+%! HTLS = freqs(BTLS, ATLS, w);
+%! [BMLS, AMLS, SigMLS] = invfreqs(H0+1e-5*Nn, w, [2 2], order, [], [], [],
[], "method", "QR");
+%! HMLS = freqs(BMLS, AMLS, w);
%! xlabel("Frequency (rad/sec)");
%! ylabel("Magnitude");
-%! plot(w,[abs(H);abs(Hh)])
+%! plot(w,[abs(H0); abs(Hh)])
%! legend('Original','Measured');
-%! err = norm(H-Hh);
+%! err = norm(H0-Hh);
%! disp(sprintf('L2 norm of frequency response error = %f',err));
--- invfreqz.m 2010-08-05 18:19:51.000000000 +0200
+++ /users/dupuis/matlab/octave-forge/signal-1.0.11/invfreqz.m 2010-10-06
08:58:41.000000000 +0200
@@ -44,9 +44,11 @@
% TODO: check invfreq.m for todo's
-function [B,A] = invfreqz(H,F,nB,nA,W,iter,tol,tr)
+function [B, A, SigN] = invfreqz(H, F, nB, nA, W, iter, tol, tr, varargin)
-if nargin < 8
+if nargin < 9
+ varargin = {};
+ if nargin < 8
tr = '';
if nargin < 7
tol = [];
@@ -57,6 +59,7 @@
end
end
end
+ end
end
@@ -60,23 +63,31 @@
end
-
% now for the real work
-[B,A] = invfreq(H,F,nB,nA,W,iter,tol,tr,'z');
+[B, A, SigN] = invfreq(H, F, nB, nA, W, iter, tol, tr, 'z', varargin{:});
+
endfunction
%!demo
-%! order = 12; % order of test filter
+%! 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
-%! [B,A] = butter(order,fc);
-%! [H,w] = freqz(B,A,n);
-%! [Bh,Ah] = invfreqz(H,w,order,order);
-%! [Hh,wh] = freqz(Bh,Ah,n);
+%! [B0, A0] = butter(order, fc);
+%! [H0, w] = freqz(B0, A0, n);
+%! Nn = (randn(size(w))+j*randn(size(w)))/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);
%! xlabel("Frequency (rad/sample)");
%! ylabel("Magnitude");
-%! plot(w,[abs(H);abs(Hh)])
+%! plot(w,[abs(H0) abs(Hh)])
%! legend('Original','Measured');
-%! err = norm(H-Hh);
+%! err = norm(H0-Hh);
%! disp(sprintf('L2 norm of frequency response error = %f',err));
------------------------------------------------------------------------------
Download new Adobe(R) Flash(R) Builder(TM) 4
The new Adobe(R) Flex(R) 4 and Flash(R) Builder(TM) 4 (formerly
Flex(R) Builder(TM)) enable the development of rich applications that run
across multiple browsers and platforms. Download your free trials today!
http://p.sf.net/sfu/adobe-dev2dev
_______________________________________________
Octave-dev mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
