FRF - PSD Calculation Using Welch's Method Producing Less Sharp Peaks for a Linear System

Hi,

I am attempting to compute the Frequency Response Function (FRF) of a linear system using Welch’s method for estimating the Power Spectral Density (PSD). However, the results I obtain lack the sharp peaks expected for a linear system, and I would like guidance on whether my code or approach is correct, or if I might be missing something fundamental.

Despite using these settings and varying the parameters, the FRF I compute does not show sharp peaks that I expect for a linear system. The peaks are relatively broad, which makes it difficult to identify resonant frequencies clearly. I am uncertain whether:

1. My choice of parameters (e.g., nfft,window_length,noverlapnfft, window\_length, noverlapnfft,window_length,noverlap) is contributing to this issue.
2. There might be a problem in my preprocessing or interpretation of the data.

Please advise, thankyou!

Fs = 8192;
window_length = 3144;
noverlap = 0.75*window_length;
nfft = 2048;

FileName = 'exp1_0';
FileStart = 1;
FileEnd = 10;
fileDirectory = 'D:\xx\';

% Define the window function
window_function = hanning(window_length);

% Initialize accumulators for averaging
sum_Pxx = 0; % Sum of auto-power spectrum of voltage
sum_Pxy = 0; % Sum of cross-power spectrum

% Loop over all test files = averaging
for filenum = FileStart:FileEnd

    % Load data from each file
    fullFileName = [fileDirectory, FileName, sprintf('%02d', filenum), '.mat'];
    fileData = load(fullFileName);

    % Extract Velocity and Voltage Data
    velocity_data = fileData.([FileName, sprintf('%02d', filenum)]).Y(4).Data; % Output: Velocity (m/s)
    voltage_data = fileData.([FileName, sprintf('%02d', filenum)]).Y(3).Data; % Input: Voltage (V)

    % Compute PSDs using Welch's method
    [Pxx, f] = pwelch(voltage_data, window_function, noverlap, nfft, Fs); % Auto-power spectrum of voltage
    [Pxy, ~] = cpsd(velocity_data, voltage_data, window_function, noverlap, nfft, Fs); % Cross-power spectrum

    % Accumulate results for averaging
    sum_Pxx = sum_Pxx + Pxx;
    sum_Pxy = sum_Pxy + Pxy;
end

% Compute FRF (H = Pxy / Pxx)
H = sum_Pxy ./ sum_Pxx;

% Convert to dB
H_mag_dB = 20 * log10(abs(H));

% Apply smoothing and moving average
window_size = 5;
smoothed_H_w = smoothdata(H_mag_dB, 'movmean', window_size);

% Plot FRF
figure;
plot(f, smoothed_H_w , 'r');
hold on;
title(['FRF of Random Test with - Window Length: ', num2str(window_length), ...
       ', Overlap: ', num2str(noverlap), ...
       ', nfft: ', num2str(nfft)]);
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
grid on;
xlim([1 300]);
hold off;

[Mathieu NOE]

hello

why not first reduce or remove the smoothing ?

and movmean will create more flatening compare to gaussian window

% Apply smoothing and moving average
window_size = 5;
smoothed_H_w = smoothdata(H_mag_dB, 'movmean', window_size);

[Mathieu NOE]

wonder if we would get the same result using the function below

(one line call and you're sone ! )

function [Txy,Cxy,f] = tfe2(x,y,nfft,Fs,window,noverlap)

%TFE Transfer Function Estimate MISO.
% x must be a vector column oriented (m x 1)
% y must be a matrix of dimensions (m x n) (n channels)
%   Txy = TFE2(X,Y,NFFT,Fs,WINDOW) estimates the transfer function of the 
%   system with input X and output Y using Welch's averaged periodogram 
%   method.  X and Y are divided into overlapping sections, each of which 
%   is detrended, then windowed by the WINDOW parameter, then zero-padded 
%   to length NFFT.  The magnitude squared of the length NFFT DFTs of the 
%   sections of X are averaged to form Pxx, the Power Spectral Density of X.
%   The products of the length NFFT DFTs of the sections of X and Y are 
%   averaged to form Pxy, the Cross Spectral Density of X and Y.  Txy
%   is the quotient of Pxy and Pxx; it has length NFFT/2+1 for NFFT even, 
%   (NFFT+1)/2 for NFFT odd, or NFFT if X or Y is complex. If you specify 
%   a scalar for WINDOW, a Hanning window of that length is used.  Fs is 
%   the sampling frequency which does not effect the transfer function 
%   estimate but is used for scaling of plots.
%
%   [Txy,F] = TFE2(X,Y,NFFT,Fs,WINDOW,NOVERLAP) returns a vector of freq-
%   uencies the same size as Txy at which the transfer function is 
%   estimated, and overlaps the sections of X and Y by NOVERLAP samples.
%
%
%   The default values for the parameters are NFFT = 256 (or LENGTH(X),
%   whichever is smaller), NOVERLAP = 0, WINDOW = HANNING(NFFT), Fs = 2, 
%   P = .95, and DFLAG = 'none'.  You can obtain a default parameter by 
%   leaving it off or inserting an empty matrix [], e.g. TFE(X,Y,[],10000).
%
%   See also PSD, CSD, COHERE
%   ETFE, SPA, and ARX in the Identification Toolbox.

% 	Author(s): M NOé, 
%       Revision: 1.0 18/04/2003

    
% 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

[p,q] = size(y);
if p~=n & q~=n
    error('x an y have incompatible dimensions');
end

if q==n         % y in wrong orientation
    y = y';
end
[p,q] = size(y);

k = fix((n-noverlap)/(nwind-noverlap));	% Number of windows
					% (k = fix(n/nwind) for noverlap=0)
index = 1:nwind;

Pxx = zeros(nfft,1);% Pxx2 = zeros(nfft,1);
Pxy = zeros(nfft,q);% Pxy2 = zeros(nfft,q);
Pyy = zeros(nfft,q);
for i=1:k
    xw = window.*x(index);
    yw = (window*ones(1,q)).*y(index,:);
    index = index + (nwind - noverlap);
    Xx = fft(xw,nfft);
    Yy = fft(yw,nfft,1);
    Xx2 = abs(Xx).^2;
    Xy2 = Yy.*(conj(Xx)*ones(1,q));
    Yy2 = abs(Yy).^2;

    Pxx = Pxx + Xx2;
    %Pxx2 = Pxx2 + abs(Xx2).^2;
    Pxy = Pxy + Xy2;
    %Pxy2 = Pxy2 + Xy2.*conj(Xy2);
    Pyy = Pyy + Yy2;

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);
    %Pxx2 = Pxx2(select);
    Pxy = Pxy(select,:);
    %Pxy2 = Pxy2(select);
    Pyy = Pyy(select,:);

else
    select = 1:nfft;
end

% set up output parameters
Txy = Pxy ./ (Pxx*ones(1,q));               % transfer function estimate 
Cxy = (abs(Pxy).^2)./((Pxx*ones(1,q)).*Pyy);% coherence function estimate 
f = (select - 1)'*Fs/nfft;                  % freq vector

end

[Mathieu NOE]

hello again

is your problem solved or do you still need our support ?

[Lucyana]

Thank you for your quick response! I tried using the function, but it didn’t work either. However, I resolved the issue by using cpsd for both Pxy and Pxx, while keeping the smooth function.

[Mathieu NOE]

if you are keen to share the data I can show you how to use the function

all the best