Clear Filters
Clear Filters

How to use PRBS signal to identify first order transfer parameters

15 views (last 30 days)
Dear All;
one of the common way to identify process model parameters is use to step test which is easy way to do but how can can use PRBS signal to do esitmate model parameters?

Answers (1)

Mathieu NOE
Mathieu NOE on 1 Sep 2023
hello
see example below
from the Bode plot you can easily extract the plant dc gain (here 10) and cut off frequency (fc = 50 Hz here, corresponding to a -3 dB gain drop vs dc gain and a phase rotation of -45 degrees)
fs = 1000; % sampling frequency (Hz)
dt = 1/fs;
%% generate the test signal
samples = 10000;
% x = rand(samples,1); % random
x = 0.5*(1+sign(randn(samples,1)));% "home made" simplified PRBS
%% plant = 1st order filter
N = 1; % filter order
fc = 50; % Hz cut off frequency
DC_gain = 10; % 1st order plant dc gain
%% generate the filtered output signal
[B,A] = butter(N,2*fc/fs);
B = B*DC_gain;
y = filter(B,A,x);
t = (0:samples-1)*dt;
figure(1)
subplot(2,1,1),plot(t(1:100),x(1:100),'b')
subplot(2,1,2),plot(t(1:100),y(1:100),'r')
%% solution 1 with tfestimate (requires Signal Processing Tbx)
% [Txy,F] = tfestimate(x,y,hanning(NFFT),NOVERLAP,NFFT,fs);
%% alternative with supplied sub function
NFFT = 512;
NOVERLAP = round(0.5*NFFT); % 0 percent overlap
[Txy,Cxy,F] = mytfe_and_coh(x,y,NFFT,fs,ones(NFFT,1),NOVERLAP);
% Txy = transfer function (complex), Cxy = coherence, F = freq vector
% Bode plots
fmax = fs/2.56; % plot only values of freq between 0 and fs/2.56 to avoid drop at Nyquist frequency
figure(2),
subplot(3,1,1),semilogx(F,20*log10(abs(Txy)));
xlim([0 fmax]);
ylabel('Mag (dB)');
subplot(3,1,2),semilogx(F,180/pi*(angle(Txy)));
xlim([0 fmax]);
ylabel('Phase (°)');
subplot(3,1,3),semilogx(F,Cxy);
xlim([0 fmax]);
ylim([0 1]);
xlabel('Frequency (Hz)');
ylabel('Coherence');
%%%%%%%%%%%%%%%%%%%%%%
function [Txy,Cxy,f] = mytfe_and_coh(x,y,nfft,Fs,window,noverlap)
% Transfer Function and Coherence Estimate
% compute PSD and CSD
window = window(:);
n = length(x); % Number of data points
nwind = length(window); % length of window
if n < nwind % zero-pad x , y if length is less than the window length
x(nwind)=0;
y(nwind)=0;
n=nwind;
end
x = x(:); % Make sure x is a column vector
y = y(:); % Make sure y is a column vector
k = fix((n-noverlap)/(nwind-noverlap)); % Number of windows
% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;
Pxx = zeros(nfft,1);
Pyy = zeros(nfft,1);
Pxy = zeros(nfft,1);
for i=1:k
xw = window.*x(index);
yw = window.*y(index);
index = index + (nwind - noverlap);
Xx = fft(xw,nfft);
Yy = fft(yw,nfft);
Xx2 = abs(Xx).^2;
Yy2 = abs(Yy).^2;
Xy2 = Yy.*conj(Xx);
Pxx = Pxx + Xx2;
Pyy = Pyy + Yy2;
Pxy = Pxy + Xy2;
end
% Select first half
if ~any(any(imag([x y])~=0)) % if x and y are not complex
if rem(nfft,2) % nfft odd
select = [1:(nfft+1)/2];
else
select = [1:nfft/2+1]; % include DC AND Nyquist
end
Pxx = Pxx(select);
Pyy = Pyy(select);
Pxy = Pxy(select);
else
select = 1:nfft;
end
Txy = Pxy ./ Pxx; % transfer function estimate
Cxy = (abs(Pxy).^2)./(Pxx.*Pyy); % coherence function estimate
f = (select - 1)'*Fs/nfft;
end

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!