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

1 view (last 30 days)
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 Comments
Smriti Mohan Shetty
Smriti Mohan Shetty on 25 Oct 2022
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!

Sign in to comment.

Accepted Answer

Mathieu NOE
Mathieu NOE on 26 Oct 2022
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 Comments

Sign in to comment.

More Answers (0)

Products


Release

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!