MSI Exercise 10 Solution

Dateiinhalt

clear all
% close all
clc

n_periods = 3;

N = 2880;
period = 1:N;


%% Multi-sine random phases

% the DFT of q_h is initialized such that it contains N zeros
Q_h = zeros(N, 1);
% the amplitude of the frequency components are set to 1 and the phases
% are set to a random value
Q_h(2:N/2) = exp(1i*2*pi*rand(N/2-1,1));
% Q_h is made conjugate symmetric
Q_h(end:-1:N/2+2) = conj(Q_h(2:N/2));
% q_h is defined using the inverse FFT
q_h = ifft(Q_h);
% q_h is rescaled in the interval [0, 6000]
q_h = 3000*(1+q_h/max(abs(q_h)));
% Q_h is updated
Q_h = fft(q_h);

%% Apply signal to system

R1a = 5e-3;
R2a = 10e-3;
R12 = 10e-3;
C1 = 1e6;
C2 = 7e5;

A = [-1/(R12*C1)-1/(R1a*C1), 1/(R12*C1); 1/(R12*C2), -1/(R12*C2)-1/(R2a*C2)];
B = [1/C1, 0; 0, 1/C2];
C = [0, 1];

t_s = 60;

sys = ss(A,B,C,[]);

qh = repmat(q_h,n_periods,1);
qd = wgn(length(qh),1,25);

figure
plot([qh,qd]);

[y,t,x] = lsim(sys,[qh,qd],(0:(N*n_periods-1))'*t_s);

figure
plot(t,x);

y_diff = abs(y(1:(n_periods-1)*N)-y(N+1:end));

figure
semilogy(y_diff);

%% Estimate transfer function

Y = fft(y((n_periods-1)*N+1:end));
U = fft(qh((n_periods-1)*N+1:end));
G = Y(2:N/2-1)./U(2:N/2-1);

figure(2);clf
semilogx(linspace(1/(2*pi*172800),1/(2*pi*240),N/2-2),20*log10(abs(G)))