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

Question

venkat ta 2023-1-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1891600-matlab-fitting-method-to-optimize-the-snr-in-the-frequency-response-curve-to-identify-high-error-fre

评论： Mathieu NOE 2023-1-12

Noise.mat

[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);

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2023-1-11

在 MATLAB Online 中打开

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);

4 个评论
显示 2更早的评论隐藏 2更早的评论

Mathieu NOE 2023-1-11

在 MATLAB Online 中打开

beam_experiment.mat

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