How to use PRBS signal to identify first order transfer parameters

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?

Answer 1

[Mathieu NOE]

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]

Problem solved ?

[Mathieu NOE]

hello again

do you mind accepting my answer (if it has fullfiled your expectations ) ? tx