# How to use PRBS signal to identify first order transfer parameters

15 views (last 30 days)
Mustafa Al-Nasser on 30 Aug 2023
Commented: Mathieu NOE on 11 Dec 2023
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?

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
Mathieu NOE on 13 Sep 2023
Problem solved ?
Mathieu NOE on 11 Dec 2023
hello again
do you mind accepting my answer (if it has fullfiled your expectations ) ? tx

### Categories

Find more on Spectral Estimation in Help Center and File Exchange

R2022a

### Community Treasure Hunt

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

Start Hunting!