![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1470886/image.png)
How to use PRBS signal to identify first order transfer parameters
10 views (last 30 days)
Show older comments
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?
0 Comments
Answers (1)
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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1470886/image.png)
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
2 Comments
Mathieu NOE
on 11 Dec 2023
hello again
do you mind accepting my answer (if it has fullfiled your expectations ) ? tx
See Also
Categories
Find more on Spectral Estimation in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!