Matlab code for data-driven frequency response estimation
The code presented below implements the method for data-driven frequency response estimation of
and reproduces the simulation results in the paper. An archive file of all functions and scripts can be downloaded from here.
Main function dd_frest
function Hh = dd_frest(ud, yd, z, n) T = n + 1; t = (1:T)'; m = size(ud, 2); p = size(yd, 2); %% preprocessing by low-rank approximation H = [moshank(ud, T); moshank(yd, T)]; [U, ~, ~] = svd(H); P = U(:, 1:m * T + n); for k = 1:length(z) %% form and solve the system of equations A = [[zeros(m * T, p); -kron(z(k) .^ t, eye(p))] P]; hg = A \ [kron(z(k) .^ t, eye(m)); zeros(p * T, m)]; Hh(:, :, k) = hg(1:p, :); end
dd_frest uses the auxiliary functions for mosaic-Hankel matrix construction
function H = moshank(w, i) if iscell(w) N = length(w); H = []; for k = 1:N, H = [H blkhank(w{k}, i)]; end else H = blkhank(w, i); end
and block-Hankel matrix construction
function H = blkhank(w, i, j) if length(size(w)) == 3 [q, N, T] = size(w); if nargin < 3 | isempty(j), j = T - i + 1; end if j <= 0, error('Not enough data.'), end H = zeros(i * q, j * N); for ii = 1:i for jj = 1:j H(((ii - 1) * q + 1):(ii * q), ... ((jj - 1) * N + 1):(jj * N)) = w(: ,:, ii + jj - 1); end end else [q, T] = size(w); if T < q, w = w'; [q, T] = size(w); end if nargin < 3 | isempty(j), j = T - i + 1; end if j <= 0, error('Not enough data.'), end H = zeros(i * q, j); for ii = 1:i H(((ii - 1) * q + 1):(ii * q), :) = w(:, ii:(ii + j - 1)); end end
Example of using dd_frest
%% simulation setup m = 2; p = 3; n = 10; Td = 200; omega = pi / 4; % frequencies over which to estimate the frequency response B = drss(n, p, m); % data-generating system h0 = freqresp(B, omega); % true frequency response e = @(hh) 100 * norm(h0 - hh) / norm(h0); % relative estimation error %% true and noisy data ud0 = rand(Td, m); yd0 = lsim(B, ud0); ud = ud0 + 0.1 * randn(size(ud0)); yd = yd0 + 0.1 * randn(size(yd0)); %% frequency response estimation with the proposed method hh_dd = dd_frest(ud0, yd0, exp(i * omega), n); e_dd = e(hh_dd)
Numerical examples in the paper
The numerical examples in this section compare the data-driven method dd_frest with the maximum-likelihood estimator implemented in the function ident, and an alternative direct data-driven estimator spa__ with the Welch method to calculate spectral densities.
function G = spa__(yd0, ud0, omega) omega = [0 omega(:)']; phi_yu = cpsd(yd0, ud0, [], [], omega); phi_uu = pwelch(ud0, [], [], omega); G = phi_yu ./ phi_uu; G = G(2:end); end
The simulation is in the errors-in-variables setup. First, we show the estimated frequency response function by the methods for a single noise realization. Then, we show the average performance of the methods obtained by Monte Carlo simulation.
Estimation of the frequency response function dd_frest_test1.m
clear all %% benchmark-system Q = [0 0 0 0.28261 0.50666]; P = [1 -1.41833 1.58939 -1.31608 0.88642]; B = ss(tf(Q, P, -1)); n = order(B); %% exact frequency response [h0, Omega] = freqresp(B); %% exact and noisy data trajectories Td = 1000; ud0 = rand(Td, 1); yd0 = lsim(B, ud0); wd0 = [ud0 yd0]; s = 0.1; wt = randn(Td, 2); wd = wd0 + s * norm(wd0) * wt / norm(wt); ud = wd(:, 1); yd = wd(:, 2); %% do the estimation with the didfferent methods hh_dd = dd_frest(ud, yd, exp(i * Omega), n); Bh = ident(wd, 1, n); hh_ident = freqresp(Bh, Omega); Bh = n4sid(iddata(yd, ud), n); hh_n4sid = freqresp(Bh, Omega); hh_spa = spa__(yd, ud, Omega); %% plot the results figure(1) plot(Omega, abs(h0(:)), '-k'), hold on, plot(Omega, abs(hh_dd(:)), '--b') plot(Omega, abs(hh_ident(:)), '-.b') plot(Omega, abs(hh_spa(:)), ':r') box off; h = gca; ax = axis; h.XTick = [ax(1) ax(2)]; h.YTick = [ax(3) ax(4)]; xlabel('frequency'), ylabel('amplitude') legend('exact', 'proposed', 'ident', 'spa', 'location', 'best', 'box', 'off') figure(2) plot(Omega, angle(h0(:)), '-k'), hold on, plot(Omega, angle(hh_dd(:)), '--b') plot(Omega, angle(hh_ident(:)), '-.b') plot(Omega, angle(hh_spa(:)), ':r') box off; h = gca; ax = axis; h.XTick = [ax(1) ax(2)]; h.YTick = [ax(3) ax(4)]; xlabel('frequency'), ylabel('phase') legend('exact', 'proposed', 'ident', 'spa', 'location', 'best', 'box', 'off')
Monte Carlo simulation dd_frest_test2.m
We show the performance of the methods as a function of the noise variance as well as the number of data points.
clear all %% benchmark-system Q = [0 0 0 0.28261 0.50666]; P = [1 -1.41833 1.58939 -1.31608 0.88642]; B = ss(tf(Q, P, -1)); n = order(B); %% Simulation parameters omega = pi / 4; h0 = freqresp(B, omega); ea = @(hh) 100 * abs(abs(h0) - abs(hh)) / abs(h0); ep = @(hh) 100 * abs(angle(h0) - angle(hh)) / angle(h0); %% True data Td = 500; ud0 = rand(Td, 1); yd0 = lsim(B, ud0); wd0 = [ud0 yd0]; %% MC simulation over different noise levels K = 11; S = linspace(0, 0.1, K); N = 100; for k = 1:K s = S(k) for i = 1:N wt = randn(Td, 2); wd = wd0 + s * norm(wd0) * wt / norm(wt); hh_dd = dd_frest(wd(:, 1), wd(:, 2), exp(j * omega), n); Bh = ident(wd, 1, n); hh_ident = freqresp(Bh, omega); hh_spa = spa__(wd(:, 2), wd(:, 1), omega); Ea_dd(i, k) = ea(hh_dd); Ep_dd(i, k) = ep(hh_dd); Ea_ident(i, k) = ea(hh_ident); Ep_ident(i, k) = ep(hh_ident); Ea_spa(i, k) = ea(hh_spa); Ep_spa(i, k) = ep(hh_spa); end end ea_dd = mean(Ea_dd); ep_dd = mean(Ep_dd); ea_ident = mean(Ea_ident); ep_ident = mean(Ep_ident); ea_spa = mean(Ea_spa); ep_spa = mean(Ep_spa); figure(1), hold on, S = S * 100; plot(S, ea_dd, 'b--') plot(S, ea_ident, 'b-.') plot(S, ea_spa, 'r:') box off; h = gca; axis([S(1) S(end) 0 ceil(max(ea_spa))]) h.XTick = [S(1) S(end)]; h.YTick = [0 ceil(max(ea_spa))]; xlabel('noise level, \%'), ylabel('error $e_a$, \%') legend('proposed', 'ident', 'spa', 'location', 'best', 'box', 'off') figure(2), hold on plot(S, ep_dd, 'b--') plot(S, ep_ident, 'b-.') plot(S, ep_spa, 'r:') box off; h = gca; axis([S(1) S(end) 0 ceil(max(ep_spa))]) h.XTick = [S(1) S(end)]; h.YTick = [0 ceil(max(ep_spa))]; xlabel('noise level, \%'), ylabel('error $e_p$, \%') legend('proposed', 'ident', 'spa', 'location', 'best', 'box', 'off') %% MC simulation over different data lengths Tmax = 1000; ud0 = rand(Tmax, 1); yd0 = lsim(B, ud0); wd0 = [ud0 yd0]; K = 10; TT = round(linspace(100, Tmax, K)); N = 100; s = 0.05; for k = 1:K T = TT(k) for i = 1:N wt = randn(T, 2); wd = wd0(1:T, :) + s * norm(wd0(1:T, :)) * wt / norm(wt); hh_dd = dd_frest(wd(:, 1), wd(:, 2), exp(j * omega), n); Bh = ident(wd, 1, n); hh_ident = freqresp(Bh, omega); hh_spa = spa__(wd(:, 2), wd(:, 1), omega); Ea_dd(i, k) = ea(hh_dd); Ep_dd(i, k) = ep(hh_dd); Ea_ident(i, k) = ea(hh_ident); Ep_ident(i, k) = ep(hh_ident); Ea_spa(i, k) = ea(hh_spa); Ep_spa(i, k) = ep(hh_spa); end end ea_dd = mean(Ea_dd); ep_dd = mean(Ep_dd); ea_ident = mean(Ea_ident); ep_ident = mean(Ep_ident); ea_spa = mean(Ea_spa); ep_spa = mean(Ep_spa); figure(3), hold on plot(TT, ea_dd(1:K), 'b--') plot(TT, ea_ident(1:K), 'b-.') plot(TT, ea_spa(1:K), 'r:') box off; h = gca; axis([TT(1) TT(end) 0 ceil(max(ea_spa))]) h.XTick = [TT(1) TT(end)]; h.YTick = [0 ceil(max(ea_spa))]; xlabel('number of samples $T$'), ylabel('error $e_a$, \%') legend('proposed', 'ident', 'spa', 'location', 'best', 'box', 'off') figure(4), hold on plot(TT, ep_dd(1:K), 'b--') plot(TT, ep_ident(1:K), 'b-.') plot(TT, ep_spa(1:K), 'r:') box off; h = gca; axis([TT(1) TT(end) 0 ceil(max(ep_spa))]) h.XTick = [TT(1) TT(end)]; h.YTick = [0 ceil(max(ep_spa))]; xlabel('number of samples $T$'), ylabel('error $e_p$, \%') legend('proposed', 'ident', 'spa', 'location', 'best', 'box', 'off')