Finding peaks and valleys, and associated indexes, of time-shifted noisy sinosidal waves

Question

Oshani 2025-9-23，20:07

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2180149-finding-peaks-and-valleys-and-associated-indexes-of-time-shifted-noisy-sinosidal-waves

评论： Mathieu NOE 2025-9-25，6:42

data.mat

I want to automate a calibration process, but I'm really struggling with the code.

I have a signal that I give my instrument (second column of attached sample data), and then two probes which record the instrument's response, independently of each other (last two columns of attached sample data). The calibration signal is a series of sine waves of varying amplitudes and periods, see below:

I need to find the following for my calibrations:

-The signal value at each peak and trough

-Whether each probe produces sinosidal waves at each period/amplitude combination

-If so, the response value for each probe at its peaks and troughs (both as individual datapoints and as the mean value) and ideally the corresponding periods and amplitudes

The first is quite easy to do and the second can be done manually without too much trouble. Where I'm running into trouble with is the third step.

I tried the findpeaks function and it's quite noisy:

I did attempt to use the MinPeakHeight argument, use a high-pass filter, and smooth my probe data using moving averages, but none of those attempts worked very well. I also tried to use a Button Down Callback Function that would store values and indexes, but that ended up being super prone to user error in terms of the exact point to click on and not very fast. Solutions like this one don't work because the two sinusoidal response aren't quite fully phase shifted, they're time-shifted by an amount of time that's not consistent either across probes or over time.

Does anyone have any suggestions for functions to look into which might help me do what I want? I'm quite stuck.

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

请先登录，再进行评论。

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

Answer 1

Mathieu NOE 2025-9-24，11:17

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2180149-finding-peaks-and-valleys-and-associated-indexes-of-time-shifted-noisy-sinosidal-waves#answer_1570562

编辑：Mathieu NOE 2025-9-24，11:23

在 MATLAB Online 中打开

data.mat

Following my first comment above , this could be another approach

extract blocks of data (you need to find the start and stop times , here I used islocalmax to find peaks
bandpass filter your probe signals , once you have measured the test signal frequency (to adapt the low and high corner frequencies of the filter
you could easily add a simple test based on the filtered probe signals :
amplitude ratio vs test signal amplitude
if the ratio is very low => this is not really a valid output
if the ratio is above 10% (for example) , again you can try to fit a sinus and extract the parameters (done in William's answer)

data=load('data.mat');

t=data.M(:,1); x=data.M(:,2); y1=data.M(:,3); y2=data.M(:,4);

samples = numel(t);

% figure

% subplot(211), plot(t,x,'-r.',t,y1,'-g.',t,y2,'-b.');

% grid on; xlabel('Time (s)'); legend('x','y1','y2');

% xlim([min(t) max(t)])

Fs = 1./mean(diff(t));

% define start/end time points of blocks of data

tmp = abs(detrend(x,1));

tmp = tmp -min(tmp); % meake tmp min equals 0

tf = islocalmax(tmp); % option 3 with islocalmax (you can also try with find peaks)

tt = t(tf);

env = tmp(tf);

figure

plot(t,tmp,tt,env,'*')

% find indexes of distant blocks

ind = find(diff(tt)>20);

tt_ind = tt(ind);

n_blocks = numel(ind);

for k= 1:n_blocks

tstop(k,1) = tt_ind(k); % end of k th block

tbefore = tt(tt<=tstop(k,1));

if k<2

tstart(k,1) = min(tbefore); % start of 1st block

else

tstart(k,1) = min(tbefore(tbefore>tstop(k-1,1))); % start of k th block

end

% add some time before and after (several seconds)

extra_time = 4; % seconds

tstop(k,1) = min(samples,tstop(k,1) + extra_time);

tstart(k,1) = max(0,tstart(k,1) - extra_time);

% now extract portion of data (k th block)

ind = (t>=tstart(k,1) & t<=tstop(k,1) );

ti = t(ind);

xx = x(ind);

yy1 = y1(ind);

yy2 = y2(ind);

% counts (input) test signal threshold crossing to get its frequency

threshold = mean(xx)+0.02; %

[t_pos,t_neg] = find_zc(ti,xx,threshold);

t_pos = unique(t_pos);

freq_input = mean(1./diff(t_pos))

% now create a bandpass filter

if freq_input<Fs/2

a = 0.5;

flow = a*freq_input;

fhigh = min(1/a*freq_input,Fs/2.56);

else

flow = 0;

fhigh = Fs/2.56;

end

[B,A] = butter(2,2/Fs*[flow fhigh]);

yy1m = mean(yy1);

yy1f = filtfilt(B,A,detrend(yy1,1));

yy1f = yy1f + yy1m;

yy2m = mean(yy2);

yy2f = filtfilt(B,A,detrend(yy2,1));

yy2f = yy2f + yy2m;

figure

plot(ti,xx,'-r.',ti,yy1f,'-g.',ti,yy2f,'-b.');

grid on; xlabel('Time (s)'); legend('x','y1','y2');

title(['Signal Freq = ' num2str(freq_input) ' Hz'])

end

freq_input = 2.3968

freq_input = 1.1509

freq_input = 0.7577

freq_input = 0.5639

freq_input = 0.2797

freq_input = 0.2189

freq_input = 0.1388

freq_input = 2.3914

freq_input = 1.1532

freq_input = 0.7566

freq_input = 0.4470

freq_input = 0.2796

freq_input = 0.2223

freq_input = 0.1384

freq_input = 2.3931

freq_input = 1.1468

freq_input = 0.7547

freq_input = 0.5647

freq_input = 0.2785

freq_input = 0.2224

freq_input = 0.1386

freq_input = 2.3883

freq_input = 1.1490

freq_input = 0.7582

freq_input = 0.5625

freq_input = 0.2796

freq_input = 0.2154

freq_input = 0.1387

freq_input = 2.3950

freq_input = 1.1506

freq_input = 0.7537

freq_input = 0.5666

freq_input = 0.2801

freq_input = 0.2230

freq_input = 0.1377

freq_input = 2.3951

freq_input = 1.1559

freq_input = 0.7595

freq_input = 0.5640

freq_input = 0.2798

freq_input = 0.2222

freq_input = 0.1242

freq_input = 2.3889

freq_input = 1.1559

freq_input = 0.7552

freq_input = 0.5666

freq_input = 0.2797

freq_input = 0.2196

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ZxP,ZxN] = find_zc(x,y,threshold)

% put data in rows

x = x(:);

y = y(:);

% 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

ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

% negative slope "zero" crossing detection, using linear interpolation

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

ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

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

Mathieu NOE 2025-9-24，12:10

编辑：Mathieu NOE 2025-9-24，12:11

在 MATLAB Online 中打开

some upgrades : added the sinus fit using from this (excellent) Fex submission (I don't have the curve fitting toolbox !)

Fminspleas - File Exchange - MATLAB Central

you can now use the results to do your tests (fit model : y = A*sin(B*x+C)+D )

hope it helps

fyi , this is the last block of data plot :

data=load('data.mat');
t=data.M(:,1); x=data.M(:,2); y1=data.M(:,3); y2=data.M(:,4);
samples = numel(t);
% figure
% subplot(211), plot(t,x,'-r.',t,y1,'-g.',t,y2,'-b.'); 
% grid on; xlabel('Time (s)'); legend('x','y1','y2'); 
% xlim([min(t) max(t)])
Fs = 1./mean(diff(t));
% define start/end time points of blocks of data 
tmp = abs(detrend(x,1));
tmp = tmp -min(tmp); % meake tmp min equals 0
tf = islocalmax(tmp); % option 3 with islocalmax (you can also try with find peaks)
tt = t(tf);
env = tmp(tf);
figure
plot(t,tmp,tt,env,'*')
% find indexes of distant blocks
ind = find(diff(tt)>20);
tt_ind = tt(ind);
n_blocks = numel(ind);
for k= 1:n_blocks
    tstop(k,1) = tt_ind(k); % end of k th block 
    tbefore = tt(tt<=tstop(k,1));
    if k<2
        tstart(k,1) = min(tbefore);  % start of 1st block 
    else
        tstart(k,1) = min(tbefore(tbefore>tstop(k-1,1))); % start of k th block
    end
    % add some time before and after (several seconds)
    extra_time = 4; % seconds
    tstop(k,1) = min(samples,tstop(k,1) + extra_time);
    tstart(k,1) = max(0,tstart(k,1) - extra_time);
    
    % now extract portion of data (k th block)
    ind = (t>=tstart(k,1) & t<=tstop(k,1) );
    ti = t(ind);
    xx = x(ind);
    yy1 = y1(ind);
    yy2 = y2(ind);
    % counts (input) test signal threshold crossing to get its frequency
    threshold = mean(xx)+0.02; % 
    [t_pos,t_neg] = find_zc(ti,xx,threshold);
    t_pos = unique(t_pos);
    freq_input = mean(1./diff(t_pos))
    
    % now create a bandpass filter
    if freq_input<Fs/2
    a = 0.5;
    flow = a*freq_input;
    fhigh = min(1/a*freq_input,Fs/2.56);
    else
    flow = 0;
    fhigh = Fs/2.56;
    end
    [B,A] = butter(2,2/Fs*[flow fhigh]);
    yy1m = mean(yy1);
    yy1f = filtfilt(B,A,detrend(yy1,1)); 
    yy1f = yy1f + yy1m;
    % sinus fit for first signal
    % use a shorter time range for the sinus fit 
    ind = (ti>=tstart(k,1)+extra_time & ti<=tstop(k,1)-extra_time );
    tfit = ti(ind);
    [yfit1,A1,B1,C1,D1] = my_sinus_fit(tfit,yy1f(ind),freq_input);% fit model : y = A*sin(B*x+C)+D
    yy2m = mean(yy2);
    yy2f = filtfilt(B,A,detrend(yy2,1)); 
    yy2f = yy2f + yy2m;
    % sinus fit for 2nd signal 
    [yfit2,A2,B2,C2,D2] = my_sinus_fit(tfit,yy2f(ind),freq_input);% fit model : y = A*sin(B*x+C)+D
    figure
    plot(ti,xx,'-r',ti,yy1f,'-g',ti,yy2f,'-b',tfit,yfit1,'-*g',tfit,yfit2,'-*b'); 
    grid on; xlabel('Time (s)'); legend('x','y1','y2','y1 fit','y2 fit');
    title(['Signal Freq = ' num2str(freq_input) ' Hz'])
    % store data 
    out_data(k,:) = [freq_input A1 B1 C1 D1 A2 B2 C2 D2];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sinus fit using fminspleas %% 
function [yfit,A,B,C,D] = my_sinus_fit(t,y,freq)
%% Solution based on fminspleas
% fit model : y = A*sin(B*x+C)+D
% There is a tool on the FEX, called fminspleas . It does all of the hard work for you.
% http://www.mathworks.com/matlabcentral/fileexchange/10093-fminspleas
% All you will need to do is define the function to be used. Here they would be something like...
funslist = {1, @(coef,xdata) cos(coef*xdata), @(coef,xdata) sin(coef*xdata)};
% Fminspleas will return 4 coefficients, in two different return variables. 
% Here, Bstart is the starting value for B. It is the ONLY one you need to provide.
%%%%%%%%%%%%% main code %%%%%%%%%%%%%%%%%
yu = max(y);
yl = min(y);
% ym = mean(y);                                             % Estimate offset
ym = (yl + yu)/2;                                             % Estimate offset
yr = (yu-yl)/2;                                             % Estimate Amplitude of 'y'
yz = y-ym;
% yz = fastsmooth(yz,4,2,1);                       % smooth data to remove noise artifacts (adjust factors)
% % Find zero-crossings
% [ZxP,ZxN] = find_zc(t,yz,0);
% per = 2*abs(ZxP(1) - ZxN(1)); % Estimate period
% freq = 1/per                                 % Estimate FREQUENCY
%%%%%%%%%%
Bstart = 2*pi*freq;
[B,ILP] = fminspleas(funslist,Bstart,t,y);
D = ILP(1);
C = atan2(ILP(2),ILP(3));
A = ILP(2)/sin( C );
yfit = A*sin(B*t+C)+D;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [ZxP,ZxN] = find_zc(x,y,threshold)
    % put data in rows
    x = x(:);
    y = y(:);
    
    % 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 
    ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
    
    % negative slope "zero" crossing detection, using linear interpolation
    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 
    ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));
end

William Rose 2025-9-24，16:11

Very nice answer by @Mathieu NOE.

Mathieu NOE 2025-9-25，6:42

tx @William Rose

请先登录，再进行评论。

Answer 2

Star Strider 2025-9-24，18:02

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2180149-finding-peaks-and-valleys-and-associated-indexes-of-time-shifted-noisy-sinosidal-waves#answer_1570585

编辑：Star Strider 2025-9-24，18:18

在 MATLAB Online 中打开

data.mat

These data are not eqsy to work with, and there is so much noise that effective filtering is impossible. Also, the data were not uniformly sampled. I resampled them to a common sampling interval before proceeding with the processing. (With the exception of the nufft function, all signal processing techniques require the data to be uinformly sampled.) My code will work for non-uniformly sampled data, however I use the resampled data for convenience throughout.

There is no 'good' way to determine if they are sinusoidal. I did a regression on a sinusiod function and then did some statistics on it, however I am not convinced that they are reliable indicators of a good fit to the sinusoidal function.

Try this --

LD = load('data.mat')

LD = struct with fields:

M: [15764×4 double]

t = LD.M(:,1); % Original Time Vector

s = LD.M(:,2); % Original Signal Vector

r = LD.M(:,[3 4]); % Original Response (Snesor) Vectors

t_stats = [mean(diff(t)); std(diff(t)); std(diff(t))/mean(diff(t))] % Statistics

t_stats = 3×1

0.1518 0.0182 0.1197

Fs = 1/t_stats(1) % Approximate Sampling Frequency

Fs = 6.5867

[yr,tr] = resample([s r], t, Fs); % Resample To Constant Sampling Frequency (Sampling Interval)

tr_stats = [mean(diff(tr)); std(diff(tr)); std(diff(tr))/mean(diff(tr))] % Statistics

tr_stats = 3×1

0.1518 0.0000 0.0000

sr = yr(:,1); % Resmapled Signal Vector

rr = yr(:,[2 3]); % Resampled Response (Sensor) Vectors

% figure

% plot(tr, sr, DisplayName='Signal')

% hold on

% plot(tr, rr(:,1), DisplayName="Sensor 1")

% plot(tr, rr(:,2), DisplayName="Sensor 2")

% hold off

% grid

% legend(Location='best')

% xlabel('Time (unit)')

% ylabel('Amplitude (unit)')

[FTs1,Fv] = FFT1(sr,tr);

figure

plot(Fv, abs(FTs1)*2)

grid

xlabel('Frequency (unit)')

ylabel('Magnitude (unit)')

title('Fourier Transform of ''Signal''')

[FTs1,Fv] = FFT1(rr,tr);

figure

plot(Fv, abs(FTs1)*2)

grid

xlabel('Frequency (unit)')

ylabel('Magnitude (unit)')

title('Fourier Transform of ''Sensors''')

% sr_filt = bandpass(sr, [0.2 0.4], Fs, ImpulseResponse='iir');

[srpks,srlocs] = findpeaks(sr, MinPeakProminence=0.05);

idx1 = find(diff([1; srlocs])>200);

fliplocs = flip(srlocs);

idx12 = find(diff([numel(tr); fliplocs])<-200);

startlocs = srlocs(idx1)-2;

endlocs = [flip(fliplocs(idx12)); numel(sr)-2]+2;

figure

hp1 = plot(tr, sr, DisplayName='Signal');

hold on

hp2 = plot(tr, rr(:,1), DisplayName="Sensor 1");

hp3 = plot(tr, rr(:,2), DisplayName="Sensor 2");

hold off

grid

hxl1 = xline(tr(startlocs)-1,'--g', DisplayName='Start');

hxl2 = xline(tr(endlocs)+1,'--r', DisplayName='End');

xlabel('Time (unit)')

ylabel('Amplitude (unit)')

legend([hp1,hp2,hp3,hxl1(1),hxl2(1)], Location='best')

figure

hp1 = plot(tr, sr, DisplayName='Signal');

hold on

hp2 = plot(tr, rr(:,1), DisplayName="Sensor 1");

hp3 = plot(tr, rr(:,2), DisplayName="Sensor 2");

hold off

grid

hxl1 = xline(tr(startlocs)-1,'--g', DisplayName='Start', LineWidth=1.5);

hxl2 = xline(tr(endlocs)+1,'--r', DisplayName='End', LineWidth=1.5);

xlabel('Time (unit)')

ylabel('Amplitude (unit)')

legend([hp1,hp2,hp3,hxl1(1),hxl2(1)], Location='best')

xlim([740 890])

title('Zoomed Segment')

●

idxofst = 20;

% figure

% for k = 1:10:numel(startlocs)

% idxrng = max(1,startlocs(k)-idxofst) : min(endlocs(k)+idxofst,numel(tr));

% figure

% plot(tr(idxrng), sr(idxrng))

% hold on

% plot(tr(idxrng), rr(idxrng,1))

% plot(tr(idxrng), rr(idxrng,2))

% hold off

% grid

% xlabel('Time (unit)')

% ylabel('Amplitude (unit)')

% title("Epoch "+k)

% end

nsl = numel(startlocs)

nsl = 49

fitfcn = @(b,t) b(1).*sin(b(2)*2*pi*t + b(3)*pi) + b(4)+b(5).*t;

f0 = 0.5;

dc = mean(sr);

figure

% % tiledlayout(numel(startlocs), size(rr,2), TileIndexing='columnmajor')

tiledlayout(5, size(rr,2), TileIndexing='columnmajor')

for k1 = 1:size(rr,2)

for k2 = 1:numel(startlocs)

% INDEX = [k1 k2]

idxrng = max(1,startlocs(k2)-idxofst) : min(endlocs(k2)+idxofst,numel(tr));

[rmin,rmax] = bounds(rr(idxrng,k1));

% Q3 = nnz(islocalmax(detrend(rr(idxrng,k1),1),MinProminence=0.01))

f0 = nnz(islocalmax(detrend(rr(idxrng,k1),2),MinProminence=0.01))/(tr(idxrng(end))-tr(idxrng(1)));

B0 = [rmax-rmin; f0; rand; dc; rand];

[B{k1,k2},nrmrsd(k1,k2)] = fminsearch(@(b)norm(fitfcn(b,tr(idxrng))-rr(idxrng,k1)), B0);

% Q4 = fitfcn(B{k1,k2},tr(idxrng))

% Q5 = rr(idxrng,k1)

% [rc(k1,k2),p(k1,k2)] = corrcoef(rr(idxrng,k1), fitfcn(B{k1,k2},tr(idxrng)));

[r,p] = corrcoef(rr(idxrng,k1), fitfcn(B{k1,k2},tr(idxrng)));

rc(k1,k2) = r(2,1);

pc(k1,k2) = p(2,1);

maxrr(k1,k2) = max(rr(idxrng,k1));

if rem(k2,9) == 0

nexttile

% figure

plot(tr(idxrng), sr(idxrng),'-b')

hold on

% plot(tr(idxrng), detrend(rr(idxrng,k1),1)+dc,'--m')

plot(tr(idxrng), rr(idxrng,k1),'-c')

plot(tr(idxrng), fitfcn(B{k1,k2},tr(idxrng)),'-r')

hold off

grid

title("Sensor "+k1+", Segment: "+k2)

end

Exiting: Maximum number of function evaluations has been exceeded - increase MaxFunEvals option. Current function value: 9.414348

sgtitle('Selected Epochs')

Bmtx = cell2mat(B);

Bmtx_1 = Bmtx(1:5,:);

Bmtx_2 = Bmtx(6:10,:);

% nrmrsd

% rc

% pc

Sensor_1 = array2table([Bmtx_1.', maxrr(1,:).', nrmrsd(1,:).', rc(1,:).', pc(1,:).'], VariableNames=["Amplitude","Frequency","Phase","Intercept","Slope","Maximum","Norm Resid","R","p"] );

Sensor_2 = array2table([Bmtx_2.', maxrr(2,:).', nrmrsd(2,:).', rc(2,:).', pc(2,:).'], VariableNames=["Amplitude","Frequency","Phase","Intercept","Slope","Maximum","Norm Resid","R","p"] );

disp(Sensor_1)

Amplitude Frequency Phase Intercept Slope Maximum Norm Resid R p __________ _________ ________ _________ ___________ _______ __________ __________ __________ 0.060917 0.32508 0.88868 16.823 0.0024979 17.03 0.32374 0.015873 0.91018 -0.0035485 0.38349 0.17768 17.09 -0.00080738 17.03 0.108 0.23928 0.061055 0.064173 0.39069 0.92592 16.759 0.0016005 17.03 0.40127 0.11088 0.3573 -0.0079663 0.46087 5.7826 17.522 -0.0026145 17.03 0.15597 0.52817 4.7749e-07 0.036817 0.29569 0.19136 17.526 -0.0021142 17.13 0.44287 0.56002 7.6555e-11 -0.074363 0.22404 0.75183 17.367 -0.0012314 17.15 0.52908 0.75477 4.0368e-26 0.022751 0.17379 4.8883 17.215 -0.00056794 17.23 2.167 0.10785 0.14175 0.060763 0.48771 0.49848 16.759 0.00059366 17.03 0.3323 0.045092 0.74851 -0.0022607 0.41524 1.1192 17.969 -0.0021475 17.05 0.1754 0.25692 0.043824 -0.0053826 0.38959 1.2111 18.432 -0.0029115 17.08 0.30828 0.35122 0.0026713 0.25989 0.49163 0.69991 16.754 0.00048321 17.15 1.6833 0.20828 0.052879 0.018124 0.47556 0.57692 18.582 -0.0026736 17.25 1.1808 0.16979 0.068434 -0.038643 0.3168 1.233 19.067 -0.0032086 17.33 2.0156 0.16713 0.053594 -0.063521 0.19097 1.1439 16.024 0.0014542 17.48 3.9599 0.15042 0.038829 0.12191 0.34345 0.56418 16.741 0.00033752 17.03 0.6256 0.14816 0.285 0.22185 0.27799 0.86523 16.737 0.0003767 17.09 1.107 0.41801 0.00072313 0.26742 0.48328 0.3763 16.752 0.00032239 17.1 1.7061 -0.065473 0.58749 0.030802 0.57676 1.9003 19.077 -0.0023641 17.28 0.81011 0.3794 0.00051899 0.35778 0.29377 0.6561 12.504 0.0049373 17.649 2.1905 0.77015 5.2803e-24 0.06271 0.27448 1.1045 16.6 0.00047003 17.7 4.4719 0.10719 0.21766 0.1222 0.19548 0.61583 16.06 0.00097536 17.841 7.3221 0.1586 0.029713 0.024314 0.16263 1.2881 16.759 0.0002176 17 0.1256 0.22802 0.10056 0.31812 0.41614 0.82789 16.743 0.00026111 17.13 1.775 0.095417 0.46069 0.021789 0.60342 1.8785 19.973 -0.0025062 17.3 0.87909 0.16281 0.17491 1.0802 0.42825 0.42292 16.758 0.00025935 17.53 6.9568 -0.0058393 0.959 0.52623 0.29434 -0.16861 25.19 -0.0064422 17.9 3.1579 0.78886 7.425e-26 -0.12236 0.26685 1.3714 33.545 -0.012548 18 6.4904 0.20164 0.018137 0.16849 0.19375 6.0906 21.719 -0.0033672 18.28 10.223 0.16193 0.026823 0.11811 0.17568 0.69096 16.785 0.00016296 17.05 0.47332 0.69463 7.9487e-09 0.045869 0.41548 0.21385 -10.035 0.018541 17.2 0.7259 0.17639 0.17025 0.035636 0.48299 8.7709 25.487 -0.0055995 17.464 1.3104 0.19452 0.10405 1.3098 0.46185 0.025964 16.749 0.00027274 17.7 7.9153 0.25615 0.021821 0.73442 0.27019 1.1606 17.004 7.2211e-05 18.3 3.6865 0.83659 1.4269e-31 0.19056 0.27098 0.033725 18.339 -0.00072946 18.42 8.3196 0.18831 0.029337 -0.077596 0.18033 0.98763 8.6239 0.0050064 18.65 13.181 0.020411 0.781 0.14262 0.3254 0.35758 16.718 0.00017864 17.097 0.68636 0.27466 0.046553 0.68744 0.89746 0.84538 16.753 0.00022174 17.446 3.9062 0.10367 0.42264 1.0913 0.6519 0.031232 16.759 0.0002202 17.68 6.3646 0.21421 0.072848 -0.071661 0.429 0.56407 14.703 0.0013296 18.08 3.7315 0.11512 0.30922 -0.12927 0.397 -3.1889 74.35 -0.029651 18.94 9.4143 0.20219 0.029507 0.29181 0.26856 27.762 83.828 -0.03364 19.19 11.697 0.23894 0.0054285 0.079197 0.20364 0.49311 19.076 -0.00093187 19.49 21.1 0.05174 0.45796 0.41061 0.17544 0.9474 16.755 0.00015765 17.3 1.5802 0.68941 1.1402e-08 -0.030323 0.89556 0.75812 19.332 -0.0010064 17.63 1.9074 0.055475 0.66847 1.6032 0.48321 0.45256 16.755 0.00022442 18.1 9.8139 0.077992 0.51796 -0.059145 0.42879 0.3463 25.325 -0.0036076 18.52 4.9365 0.07619 0.50176 0.78337 0.31722 0.53493 11.813 0.0024672 19.54 10.857 0.47497 7.1608e-08 -0.10586 0.33998 0.98629 20.492 -0.0013305 19.847 15.434 0.072391 0.40407 1.0154 0.1609 0.53973 18.048 -0.00034723 21.75 23.927 0.40415 1.0746e-09

disp(Sensor_2)

Amplitude Frequency Phase Intercept Slope Maximum Norm Resid R p __________ _________ _________ _________ ___________ _______ __________ __________ __________ 0.025282 0.17838 1.0686 16.814 0.0035997 17.1 0.098206 0.76392 2.8533e-11 -0.0088895 0.13456 2.6482 16.654 0.0040947 17.13 0.037123 0.94694 2.8698e-31 0.0066364 0.35764 1.8729 16.345 0.0057916 17.28 0.083302 0.87795 9.1797e-24 0.058683 0.23014 0.27825 16.868 0.0016399 17.2 0.39202 -0.045832 0.68644 0.21254 0.31989 0.99023 16.752 0.0015573 17.23 1.714 -0.035517 0.7063 0.23062 0.20509 0.5069 16.822 0.00085202 17.18 1.6418 0.58914 5.6369e-14 0.10011 0.14709 1.7222 15.648 0.0040681 17.25 1.0002 0.62447 1.2998e-21 0.060758 0.32517 1.1533 16.754 0.00092327 17.15 0.33156 0.042393 0.76311 0.11074 0.41583 0.1742 16.752 0.00079024 17.13 0.64836 -0.035223 0.78579 0.015093 0.13059 0.11068 18.289 -0.002696 16.98 0.090788 0.77112 3.6352e-15 0.0045602 0.63663 0.74701 19.301 -0.004256 17.08 0.2615 0.50834 4.9951e-07 0.16618 0.36766 0.0044297 16.757 0.00054429 17.15 1.3169 0.023942 0.79864 0.29528 0.27434 0.46919 16.757 0.00065102 17.33 2.4959 0.0045768 0.95814 0.68174 0.18034 0.29299 16.722 0.00054188 17.43 6.8649 0.011904 0.87086 0.15402 0.34448 0.1444 16.767 0.00031964 17.05 0.79092 0.16255 0.24025 0.1169 0.29937 0.25028 16.764 0.00021649 16.98 0.61994 0.18943 0.14032 0.15316 0.48414 1.1823 16.731 0.0002131 16.95 0.91277 0.21002 0.078757 0.19376 0.57695 0.58236 16.752 0.00019094 16.98 1.2132 0.15374 0.17333 -0.015889 0.31736 1.3931 20.497 -0.0037393 17.25 0.86872 0.03439 0.714 0.091956 0.19024 0.57791 19.057 -0.0019309 17.48 1.1645 0.55331 4.1141e-12 0.39848 0.13559 0.90534 28.34 -0.010829 17.73 1.6373 0.9227 6.3822e-79 0.03647 0.16262 0.57519 16.758 0.00031593 17.13 0.20519 0.091698 0.51373 0.11347 0.13864 0.32699 16.752 0.00026836 17.1 0.45694 0.7341 1.131e-11 0.16686 0.13041 0.29001 16.751 0.00023194 17.08 0.70326 0.75211 4.0272e-14 0.018116 0.23098 1.0207 3.0418 0.011599 17.15 0.50247 0.57647 2.1964e-08 0.069065 0.29404 1.8679 24.392 -0.0058441 17.35 1.2138 0.44592 5.2631e-07 -0.094306 0.24815 4.0103 28.516 -0.0087403 17.53 2.1156 0.40545 8.8831e-07 -0.10584 0.1948 1.6872 30.578 -0.0098351 17.9 6.0271 0.24141 0.00087317 0.059017 0.17544 0.53657 16.752 0.00025873 17.15 0.31158 -0.012832 0.92734 0.12145 0.1387 0.65328 16.744 0.00023526 17.15 0.42784 0.862 2.3378e-19 0.19489 0.24179 1.1421 16.727 0.00024594 17.15 1.1855 0.040321 0.73849 0.25969 0.32116 0.65857 16.752 0.00020249 17.15 1.7518 0.0061061 0.95713 1.0091 0.29444 1.2065 16.729 0.00023099 17.53 7.056 0.38383 2.107e-05 0.14146 0.25463 -0.70008 15.117 0.0012208 17.85 3.2982 0.30488 0.00034115 2.6706 0.14638 0.15013 16.752 0.00011713 18.13 20.572 0.75602 4.4765e-36 0.085462 0.16251 0.91227 16.756 0.00011486 17 0.33874 0.64325 2.0526e-07 0.30477 0.27717 0.12413 16.758 0.00012716 17.1 1.7253 0.067287 0.60332 0.42654 0.12076 0.50444 16.739 0.00018637 17.28 1.7996 0.6811 6.3083e-11 0.56042 0.21411 0.40357 16.752 0.00024697 17.38 3.6093 0.10358 0.36053 0.10854 0.31743 1.0971 18.342 -0.00062374 17.83 3.2536 0.25664 0.0054202 -0.051044 0.27396 0.74602 2.0966 0.0076086 18.13 5.3897 -0.16943 0.050339 -0.37561 0.17615 -0.013568 63.264 -0.02265 18.8 14.069 0.32544 1.6141e-06 0.20776 0.48803 0.12307 16.766 0.00011407 17.1 1.0896 0.1634 0.24238 0.41277 0.27765 0.59439 16.732 0.00017532 17.28 2.4387 -0.0011907 0.99267 0.59322 0.2603 0.49348 16.737 0.00017422 17.43 3.6562 0.064343 0.59397 0.031463 0.42814 0.2067 16.893 0.00010073 17.4 2.1171 0.1047 0.35535 -0.17194 0.31769 0.87049 25.159 -0.0035187 18.145 4.5629 0.3035 0.00092688 -0.028999 0.34242 -1.0438 27.21 -0.0043195 18.51 7.1958 0.24907 0.0035794 0.12394 0.17301 2.0674 16.928 -3.9214e-06 19.09 16.159 0.078403 0.25686

function [FTs1,Fv] = FFT1(s,t)

% One-Sided Numerical Fourier Transform

% Arguments:

% s: Signal Vector Or Matrix

% t: Associated Time Vector

t = t(:);

L = numel(t);

if size(s,2) == L

s = s.';

end

Fs = 1/mean(diff(t));

Fn = Fs/2;

NFFT = 2^nextpow2(L);

FTs = fft((s - mean(s)) .* hann(L).*ones(1,size(s,2)), NFFT)/sum(hann(L));

Fv = Fs*(0:(NFFT/2))/NFFT;

% Fv = linspace(0, 1, NFFT/2+1)*Fn;

Iv = 1:numel(Fv);

Fv = Fv(:);

FTs1 = FTs(Iv,:);

end

EDIT -- Re-ran code.

.

2 个评论
显示无隐藏无

William Rose 2025-9-24，19:43

Very nice by @Star Strider

Star Strider 2025-9-24，21:04

Thank you!

I'm not sure what Oshani wants. The sensors are not creating an accurate representation of the signal, so I have no idea what the calibration is going to do or what the actual objective is here. It may be that a relatively high-frequency sine is not appropriate for the instrumentation, and a step might be better.

I started working on this last night, however my computer crashed -- unusual for Ubuntu -- and since it ws late, I just gave up for the night. This took seeral hours to analyse.

请先登录，再进行评论。

Answer 3

William Rose 2025-9-24，8:00

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2180149-finding-peaks-and-valleys-and-associated-indexes-of-time-shifted-noisy-sinosidal-waves#answer_1570554

在 MATLAB Online 中打开

data.mat

@Oshani,

You write that you want (and I have added numbers to your three desires),

1. The signal value at each peak and trough

2. Whether each probe produces sinosidal waves at each period/amplitude combination

3. If so, the response value for each probe at its peaks and troughs (both as individual datapoints and as the mean value) and ideally the corresponding periods and amplitudes

And you write "The first is quite easy to do and the second can be done manually without too much trouble. Where I'm running into trouble with is the third step."

I am glad that you can determine manually and without much trouble whether the waves are sinusoidal. That is beyond me. I would have to use function thd() or something like it to quantify sinusoidal-ness.

The interval between samples varies, which will complicate the analysis. The plots below show the three signals (top) and the sampling interval (middle). The input signal is labelled "x" and the outputs are "y1" and "y2". I will discuss the bottom panel later.

data=load('data.mat');
t=data.M(:,1); x=data.M(:,2); y1=data.M(:,3); y2=data.M(:,4);
figure
subplot(311), plot(t,x,'-r.',t,y1,'-g.',t,y2,'-b.'); 
grid on; xlabel('Time (s)'); legend('x','y1','y2'); 
subplot(312), plot(t(2:end),diff(t),'-k.'); 
grid on; xlabel('Time (s)'); ylabel('Sampling Interval (s)')

If the long intervals only occur between the pulses, then it would not be a big problem, but that is not the case. For example, if you zoom in on the time period 2020<=t<=2060, you will see that the long time intervals sometime occur in the middle of a sinusoidal pulse.

The bottom panel is a zoomed-in view of the first pulse. The plot shows that the signal is not sampled rapidly enough to determine if it is sinusoidal or not.

subplot(313)

plot(t,x,'-r.',t,y1,'-g.',t,y2,'-b.');

grid on; xlabel('Time'); legend('x','y1','y2');

xlim([73,76]); ylim([16.6,17.4])

The bottom panel shows 10 data points for 4 cycles of oscillation, i.e. the frerquency is 0.4 cycles per sample. If this wave is non-sinusoidal, i.e. if it has power at harmonics of the fundamental frequency, then those harmonics will be above the Nyquist frequency (0.5 cycles/sample), so they will be aliased down to lower frequencies, in the range 0 to 0.5 cycles/sample, and you will not be able to un-alias them or do a reliable analysis of the signal.

Can you resample the signals faster and at a constant rate?

To accomplish your third aim, which is to find the period and amplitude of the signal and probe responses for the sinusoidal pulses, I recommend fitting a sinusoid to the signal pulse and to the pulses recorded by the probes. There is no doubt a clever way to determine where the pulses start and end*, but for the sake of simlicity I will simply determine the pulse start and end times by inspecting the plot. Then I wil fit the data in those windows. Start and end times of first four pulses and last four pulses:

% Enter pulse start, end times, determined by inspecting the plots
Ts=[73.8,112,152,193.9,2213.8,2257.3,2308.2,2362.7];  
% pulse start times (<= time of first point in each pulse)
Te=[75.6,115.6,157.4,201.1,2221,2271.7,2326.4,2391.4];  
% pulse end times (>= time of last point in each pulse)

The next chunk of code is a funciton that fits a sinusoid to x,y data:

function [f,gof] = fitSinusoid(x,y)
% FITSINUSOID Fit sinusoid to data
%   Model: y=a*sin(b*x+c)+d
%   Make a good initial guess for parameters in order to increase
%   chance of getting a good fit.
a0=(max(y)-min(y))/2;           % amplitude initial guess
nzc=length(find(diff(sign(y-mean(y)))==-2)); 
% number of negative-going zero crossings
b0=2*pi*nzc/(max(x)-min(x));    % radian freq initial guess
c0=0;                           % phase initial guess
d0=mean(y);                     % offset initial guess
myModel=@(a,b,c,d,x) a*sin(b*x+c)+d;
[f1,gof1]=fit(x,y,myModel,...
    'StartPoint', [a0,b0,c0,d0], ...
    'Lower', [0, 0, 0, -Inf], ...
    'Upper', [Inf, Inf, max(x)-min(x), Inf]);
[f2,gof2]=fit(x,y,myModel,...
    'StartPoint', [2*a0,b0,c0,d0], ...
    'Lower', [0, 0, 0, -Inf], ...
    'Upper', [Inf, Inf, max(x)-min(x), Inf]);
if gof1.sse<=gof2.sse
    f=f1; gof=gof1;
else
    f=f2; gof=gof2;
end
end

You may be wondering why I didn't just use the built-in fittype 'sin1'. The reason is that 'sin1' does not include a mean offset, as far as I can tell.

The next chunk of code fits the signal and probe repsonses during each pulse, and displays the best-fit amplitude and frequency and r^2 on the console. It also displays plots of the data and the fits.

% Analyze each pulse

Np=length(Ts); % number of pulses

for i=1:Np

[fx,gofx] =fitSinusoid(t(t>=Ts(i) & t<=Te(i)),x(t>=Ts(i) & t<=Te(i)));

[fy1,gofy1]=fitSinusoid(t(t>=Ts(i) & t<=Te(i)),y1(t>=Ts(i) & t<=Te(i)));

[fy2,gofy2]=fitSinusoid(t(t>=Ts(i) & t<=Te(i)),y2(t>=Ts(i) & t<=Te(i)));

fprintf('Pulse %d:\n',i);

fprintf(' Signal: ampl=%.3f, freq=%.3f, r^2=%.2f.\n',...

fx.a,fx.b/(2*pi),gofx.rsquare)

fprintf(' Probe1: ampl=%.3f, freq=%.3f, r^2=%.2f.\n',...

fy1.a,fy1.b/(2*pi),gofy1.rsquare)

fprintf(' Probe2: ampl=%.3f, freq=%.3f, r^2=%.2f.\n',...

fy2.a,fy2.b/(2*pi),gofy2.rsquare)

figure

subplot(311)

plot(fx,t(t>=Ts(i) & t<=Te(i)),x(t>=Ts(i) & t<=Te(i)))

xlabel('Time'); grid on

subplot(312)

plot(fy1,t(t>=Ts(i) & t<=Te(i)),y1(t>=Ts(i) & t<=Te(i)))

xlabel('Time'); grid on

subplot(313)

plot(fy2,t(t>=Ts(i) & t<=Te(i)),y2(t>=Ts(i) & t<=Te(i)))

xlabel('Time'); grid on

end

Pulse 1:

Signal: ampl=0.252, freq=2.407, r^2=1.00.

Probe1: ampl=0.000, freq=0.000, r^2=NaN.

Probe2: ampl=0.000, freq=0.000, r^2=0.00.

Pulse 2:

Signal: ampl=0.249, freq=1.151, r^2=1.00.

Probe1: ampl=0.016, freq=0.291, r^2=0.66.

Probe2: ampl=0.030, freq=0.000, r^2=0.00.

Pulse 3:

Signal: ampl=0.250, freq=0.756, r^2=1.00.

Probe1: ampl=0.012, freq=0.379, r^2=0.23.

Probe2: ampl=0.050, freq=0.000, r^2=0.00.

Pulse 4:

Signal: ampl=0.250, freq=0.562, r^2=1.00.

Probe1: ampl=0.003, freq=0.421, r^2=0.01.

Probe2: ampl=0.002, freq=0.278, r^2=0.03.

Pulse 5:

Signal: ampl=3.999, freq=0.567, r^2=1.00.

Probe1: ampl=0.125, freq=0.424, r^2=0.03.

Probe2: ampl=0.516, freq=0.000, r^2=0.00.

Pulse 6:

Signal: ampl=4.000, freq=0.280, r^2=1.00.

Probe1: ampl=1.551, freq=0.279, r^2=0.86.

Probe2: ampl=0.875, freq=0.140, r^2=-1.08.

Pulse 7:

Signal: ampl=3.995, freq=0.220, r^2=0.99.

Probe1: ampl=1.859, freq=0.219, r^2=0.87.

Probe2: ampl=0.164, freq=0.220, r^2=0.22.

Pulse 8:

Signal: ampl=0.282, freq=0.139, r^2=0.14.

Probe1: ampl=2.449, freq=0.138, r^2=0.93.

Probe2: ampl=1.385, freq=0.138, r^2=0.76.

The plots show that the fits are not always the best, but maybe with some tweaking, you can get this approach to be useful. Maybe use fmincon() instead of fit(). Also, the script reports r^2=-1.08 for pulse 6. That shouldn't happen, so it merits further investigation.

Good luck with your analysis.

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

Mathieu NOE 2025-9-24，9:01

在 MATLAB Online 中打开

data.mat

dear all

IMHO the test signals have a too high frequency , the probe will not keep with the firsts sine bursts that are obviously beyond their frequency bandwidth

maybe you should redo the test with lower freq bursts (and that would also increase the number of samples per period )

to test the sinus like shape maybe we could also use first some bandpass filtering of the probe signals , based on the frequency measured on the input (test) signal. frequency can be measured from the "zero" crossing rate of the signal (using the mean of the signal as the threshold value) - see second subplot here

data=load('data.mat');

t=data.M(:,1); x=data.M(:,2); y1=data.M(:,3); y2=data.M(:,4);

figure

subplot(211), plot(t,x,'-r.',t,y1,'-g.',t,y2,'-b.');

grid on; xlabel('Time (s)'); legend('x','y1','y2');

xlim([min(t) max(t)])

% counts signal crossing

threshold = mean(x)+0.1; %

[t0_pos,t0_neg] = find_zc(t,x,threshold);

freq_input = [0;1./diff(t0_pos)];

% discard values below 0.05 Hz (generated by distant t0_pos belonging to separate burst)

ind = (freq_input<0.05);

freq_input(ind) = [];

t0_pos(ind) = [];

subplot(212), plot(t0_pos,freq_input,'*');

grid on; xlabel('Time (s)'); ylabel('Signal freq (Hz)')

xlim([min(t) max(t)])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [ZxP,ZxN] = find_zc(x,y,threshold)

% put data in rows

x = x(:);

y = y(:);

% 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

ZxP = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

% negative slope "zero" crossing detection, using linear interpolation

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

ZxN = ZeroX(x(ix),y(ix),x(ix+1),y(ix+1));

end

请先登录，再进行评论。

Finding peaks and valleys, and associated indexes, of time-shifted noisy sinosidal waves

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

回答（3 个）

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

2 个评论
显示无隐藏无

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Finding peaks and valleys, and associated indexes, of time-shifted noisy sinosidal waves

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

回答（3 个）

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

2 个评论 显示 无隐藏 无

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

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

2 个评论
显示无隐藏无

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