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

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!

[Mathieu NOE]

hello

what is the issue ? the some seems fine at first look

please attach the mat file (paperclip button) if you want someone to test your code

all the best

[Sai]

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]

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]

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

[Smriti Mohan Shetty]

Hi, thanks for this! I do also want to use it to find amplitude, could i do that?

[Mathieu NOE]

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]

Hello,

I keep getting an error for the second line in the code "data1=data1.theta"(that it is an unrecognized variable. I tried changing it back to "data1=data1.pos" as well, that doesn't work either. Could you attach the mat files so I can see if I'm missing something here.

[Mathieu NOE]

hello

this is the mat file - as you posted above !

the code still works fine with this file and with "data1=data1.theta"

[Smriti Mohan Shetty]

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]

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