Matlab fitting method to optimize the SNR in the frequency response curve to identify high error frequencies

[stmfile,stmpath]=uigetfile('*mat','pick the mat file');

File       = fullfile(char(stmpath),char(stmfile));
load(File);   

[~,fileName,~] = fileparts(char(File)); % e.g., file is 'dp600_2_layers_L50.xlsx'



semilogx(No_smooth_x,'r-');hold on

grid on
grid minor

xlabel('Frequency(Hz)','FontSize',20)

set(gca,'FontSize',20);

Answer 1

[Mathieu NOE]

hello

this will reduce your plot noise but maybe you should improve the measurement method first ?

load('Noise.mat');

smooth_x = smoothdata(No_smooth_x,'gaussian',500);
% keep original data below f = 500 Hz
smooth_x(1:500) = No_smooth_x(1:500);

semilogx(No_smooth_x,'r-');hold on
semilogx(smooth_x,'b-');hold on
grid on
grid minor
xlabel('Frequency(Hz)','FontSize',20)
set(gca,'FontSize',20);

[venkat ta]

Yes. The id is to optimize the measured data SNR at high frequency and resonance frequency by comparing the fitted curve (exactly as you smoothed it) and the noisy curve then to increase the power of the noise where error > 3dB.

seems smoothdata is also considered like fitted curve

[Mathieu NOE]

When I do ID with time data I like also to plot the coherence so I know if my data is good or not in agiven frequency range. It can be simply to optimize the signal strength or increase the record duration to have more averages

here some demo code :

time data in mat file attached

data = load('beam_experiment.mat');
x = transpose(data.x); %input
y = transpose(data.y); %output
fs = data.fs;   % sampling frequency

NFFT = 2048;
NOVERLAP = round(0.75*NFFT);  % 75 percent overlap

%% solution 1 with tfestimate (requires Signal Processing Tbx)
% [Txy,F] = tfestimate(x,y,hanning(NFFT),NOVERLAP,NFFT,fs);

%% alternative with supplied sub function 
[Txy,Cxy,F] = mytfe_and_coh(x,y,NFFT,fs,hanning(NFFT),NOVERLAP); 
% Txy = transfer function (complex), Cxy = coherence, F = freq vector

% Bode plots
figure(1),
subplot(3,1,1),plot(F,20*log10(abs(Txy)));
ylabel('Mag (dB)');
subplot(3,1,2),plot(F,180/pi*(angle(Txy)));
ylabel('Phase (°)');
subplot(3,1,3),plot(F,Cxy);
xlabel('Frequency (Hz)');
ylabel('Coh');

%%%%%%%%%%%%%%%%%%%%%%%
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

[venkat ta]

Nice example to follow, thx

[Mathieu NOE]

My pleasure !