Help with finding frequency, amplitude and error of a signal that I'm uploading from a simulink model

Question

Smriti Mohan Shetty 2022-10-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1830313-help-with-finding-frequency-amplitude-and-error-of-a-signal-that-i-m-uploading-from-a-simulink-mode

评论： Smriti Mohan Shetty 2022-11-3

data1 = load("pos.mat");
data1 = data1.pos;
fs = 100;
x = data1(1,:);
y = data1(2,:);
t = 1.5;
lnsig=length(y);
midline = mean(y)*ones(1,lnsig);
figure(1)
plot(x,y,'b',x,t*ones(1,lnsig),'g',x,midline,'r')
tite('target position: 1.5 radians')
xlabel('Time[s]')
ylabel('Position[radians]')
legend('y vs x','target','midline')
error = y(end)-t;
[pxx,f] = pwelch(y,fs)
fhz = f/pi*fs*2;
[mxp idx] = max(pxx);
frequency = fhz(idx);
period = 1/frequency;
%%%%%%%%%%%%%%%%%%%%%%%%%%

Let me know if there are any other codes that would work. Thanks!

3 个评论
显示 1更早的评论隐藏 1更早的评论

Sai 2022-10-21

Can you provide me with the "pos.mat" file so that I can try it from my end and the detailed description of the issue

Smriti Mohan Shetty 2022-10-25

Thank you for your replies!

I have attached the code as well as the pos.mat file along with a matlab simulink model-position.mdl. I have a matlab simulink model(position.mdl) with one of the outputs going into the pos.mat file. I need the code above to be able to find the frequency,amplitude of the data in the pos.mat file each time i run the code with different values in certain blocks in the simulink model(pos reference and gain values).

Hope this helps!

请先登录，再进行评论。

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

Answer 1

Mathieu NOE 2022-10-26

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1830313-help-with-finding-frequency-amplitude-and-error-of-a-signal-that-i-m-uploading-from-a-simulink-mode#answer_1083823

在 MATLAB Online 中打开

hello again

if your data are "clean" sinewaves but the number of samples is limited then computing the period via the time diffeence between zero crossing points (or any non zro threshold as well) will give you a better precision than a fft based method (like pwelch) . Remember that the fft frequency vector resolution is given by df = Fs / nfft (where nfft is the length of the samples used in your fft computation , so here your entire signal length for best accuracy)

you can compare both values from the updated code below

pwelch method :

frequency = 1.1719 +/- 0.3906 (quite coarse ! )

period = 0.8533 +/- 0.32 (idem)

zero crossing method

period2 = 0.9225 (fine resolution : method depends on noise level and accuracy of zc linear interpolation method)

BTW, your data contained only a field theta but not pos

there was also a small mistake when you do the frequency vector computation

data1 = load("pos.mat");
% data1 = data1.pos;
data1 = data1.theta;
fs = 100;
x = data1(1,:);
y = data1(2,:);
t = 1.5;
lnsig=length(y);
midline = mean(y)*ones(1,lnsig);
error = y(end)-t;
[pxx,f] = pwelch(y,fs);
% fhz = f/pi*fs*2;
fhz = (f/pi)*(fs/2);  % correction here
[mxp, idx] = max(pxx);
frequency = fhz(idx);
period = 1/frequency
% zero crossing method -alternative-
zct = find_zc(x,y,0);
period2 = mean(diff(zct))
figure(1)
plot(x,y,'b',x,t*ones(1,lnsig),'g',x,midline,'r',zct,zeros(1,numel(zct)),'.k','markersize',50)
title('target position: 1.5 radians')
xlabel('Time[s]')
ylabel('Position[radians]')
legend('y vs x','target','midline','zc points')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function zct = find_zc(x,y,threshold)
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
zct = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end

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

Mathieu NOE 2022-10-27

在 MATLAB Online 中打开

hello again

see new code below , for both methods. I replaced the pwelch method with the simple windowed fft (didn't have the time to dig again into scaling factor for pwelch function... so I prefered to use a piece of code I'm trusting

fft method :

period1 = 1.0600

amplitude1 = 1.5029

zc method :

period2 = 0.9225

amplitude2 = 1.5702

data1 = load("pos.mat");
% data1 = data1.pos;
data1 = data1.theta;
fs = 100;
x = data1(1,:);
y = data1(2,:);
t = 1.5;
lnsig=length(y);
midline = mean(y)*ones(1,lnsig);
error = y(end)-t;
% [pxx,f] = pwelch(y,fs);
% fhz = (f/pi)*(fs/2);  % correction here
% [mxp, idx] = max(pxx);
[fhz,fft_spectrum] = do_fft(x(:),y(:));
[amplitude1, idx] = max(fft_spectrum);
frequency = fhz(idx);
period1 = 1/frequency
amplitude1
% zero crossing method -alternative for clean periodic signals-
zct = find_zc(x,y,0);
period2 = mean(diff(zct))
amplitude2 = 0.5*(max(y) - min(y))
figure(1)
plot(x,y,'b',x,t*ones(1,lnsig),'g',x,midline,'r',zct,zeros(1,numel(zct)),'.k','markersize',50)
title('target position: 1.5 radians')
xlabel('Time[s]')
ylabel('Position[radians]')
legend('y vs x','target','midline','zc points')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function zct = find_zc(x,y,threshold)
% positive slope "zero" crossing detection, using linear interpolation
y = y - threshold;
zci = @(data) find(diff(sign(data))>0);  %define function: returns indices of +ZCs
ix=zci(y);                      %find indices of + zero crossings of x
ZeroX = @(x0,y0,x1,y1) x0 - (y0.*(x0 - x1))./(y0 - y1); % Interpolated x value for Zero-Crossing 
zct = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [freq_vector,fft_spectrum] = do_fft(time,data)
dt = mean(diff(time));
Fs = 1/dt;
nfft = length(data); % maximise freq resolution => nfft equals signal length
% window : hanning
window = hanning(nfft);
cor_coef = length(window)/sum(window);
% fft scaling
% fft_spectrum = abs(fft(data))/nfft;
fft_spectrum = abs(fft(data.*window))*2*cor_coef/nfft;
% one sidded fft spectrum  % Select first half 
    if rem(nfft,2)    % nfft odd
        select = (1:(nfft+1)/2)';
    else
        select = (1:nfft/2+1)';
    end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end

Smriti Mohan Shetty 2022-11-3

EXPnew_Gain0.8_position_1.5.mat

Is there a reason for why the same code does not work for this attached file? It's just another variation of the same data.

Smriti Mohan Shetty 2022-11-3

Thanks! I was able to change it and run the file!

请先登录，再进行评论。

Help with finding frequency, amplitude and error of a signal that I'm uploading from a simulink model

3 个评论
显示 1更早的评论隐藏 1更早的评论

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Help with finding frequency, amplitude and error of a signal that I'm uploading from a simulink model

3 个评论 显示 1更早的评论隐藏 1更早的评论

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论

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