finding peak points while signal processing

clear

% Define the mass, damping, and stiffness matrices

load C K

load M.mat

% Define the input force as a sinusoidal function of frequency w in Hz

f = linspace(0, 40, 100); % linearly spaced frequency vector in Hz

w = 2*pi*f; % angular frequency vector in rad/s

F = @(t) [sin(2*pi*f(1)*t); zeros(3, 1)];

% Define the initial conditions

x0 = zeros(4, 1);

v0 = zeros(4, 1);

y0 = [x0; v0];

ind = 100;

% Define the time span

tspan = [0 20];

% Preallocate arrays to store the magnitude and phase of the response

mag = zeros(size(w));

phase = zeros(size(w));

X = zeros(length(w), 4);

% Loop over the frequency vector and compute the response at each frequency

for i = 1:length(w)

% Define the input force at the current frequency

F = @(t) [sin(w(i)*t); zeros(3, 1)];

% Solve the differential equation to obtain the steady-state response

[~, y] = ode45(@(t, y) vibration_equation(t, y, M, C, K, F), tspan, y0);

x = y(end, 1:4)';

% Compute the magnitude, phase, and displacement of the response at the current frequency

mag(i) = abs(x(1));

phase(i) = angle(x(1)) * 180/pi;

X(i,:) = x';

mpp = max(mag(1,:))/30;

[Ypk,Xpk,Wpk,Ppk] = findpeaks(mag(1,:),'MinPeakProminence',mpp, 'WidthReference','halfheight');

Xpk = Xpk + min(size(mag)) - 1;

[~,mxidx] = maxk(Ppk,4); % to get top 4 peak prominences

Ypkc{i} = Ypk(mxidx); % peak amp

Xpkc{i} = Xpk(mxidx); % peak locations

Wpkc{i} = Wpk(mxidx); % peak widths

Ppkc{i} = Ppk(mxidx); % peak prominences

freq{i} = f(Xpkc{i}); % Frequency

end

figure;

subplot(2, 1, 1);

%plot(f, mag);

plot(f,mag,f(Xpk),Ypk,'dr');

xlabel('Frequency (Hz)');

ylabel('Magnitude');

ylim([0 3.5/10000]);

title('Frequency Response Function');

grid on;

subplot(2, 1, 2);

plot(f, mod(phase+180,360)-180,f(Xpk),Ypk,'dr');

xlabel('Frequency (Hz)');

ylabel('Phase (deg)');

ylim([-180, 180]);

grid on;

Ypkc = cell(1,4);

Xpkc = cell(1,4);

Wpkc = cell(1,4);

Ppkc = cell(1,4);

for i = 1:4

% Find the minimum peak prominence as a fraction of the maximum response

mpp(1,i) = max(X(:,i))/40;

% Find the peaks in the response signal

[Ypk,Xpk,Wpk,Ppk] = findpeaks(X(:,i), 'MinPeakProminence', mpp(1,i), 'WidthReference', 'halfheight');

% Adjust the peak locations to account for the starting index of the response signal

Xpk = Xpk + 1;

% Find the top 4 peak prominences and store the peak information

[~, mxidx] = maxk(Ppk, 4);

Ypkc{i} = Ypk(mxidx);

Xpkc{i} = Xpk(mxidx);

Wpkc{i} = Wpk(mxidx);

Ppkc{i} = Ppk(mxidx);

end

% Plot the magnitude, phase, and displacement of the FRF in Hz

figure;

subplot(2, 2, 1);

plot(f, abs(X(:,1)), 'b');

hold on;

plot(f(Xpk), Ypk, 'dr');

xlabel('Frequency (Hz)');

ylabel('Displacement (m)');

ylim([0 3.5/10000]);

title('X1');

grid on;

subplot(2, 2, 2);

plot(f, abs(X(:,2)),'b');

hold on;

plot(f(Xpk), Ypk, 'dr');

xlabel('Frequency (Hz)');

ylabel('Displacement (m)');

ylim([0 3.5/10000]);

title('X2');

grid on;

subplot(2, 2, 3);

plot(f, abs(X(:,3)),'b');

hold on;

plot(f(Xpk), Ypk, 'dr');

xlabel('Frequency (Hz)');

ylabel('Displacement (m)');

ylim([0 3.5/10000]);

title('X3');

grid on;

subplot(2, 2, 4);

plot(f, abs(X(:,4)),'b');

hold on;

plot(f(Xpk), Ypk, 'dr');

xlabel('Frequency (Hz)');

ylabel('Displacement (m)');

ylim([0 3.5/10000]);

title('X4');

grid on;

clear C damp_ratio F i ind K K1_Th K2_Th M Natural_freq tspan mxidx

function dydt = vibration_equation(t, y, M, C, K, F)

% Compute the acceleration

x = y(1:4);

v = y(5:8);

a = M \ (F(t) - C*v - K*x);

% Compute the derivative of the state vector

dydt = [v; a];

end

Right now, I am getting wrong number of peaks at wrong frequancy can anyone help me with this bug.

[AL]

Good morning @Mathieu NOE,

The modification which you have suggested are working perfectly and I really appricate youe efforts. Have a wonderful day.

-AL.

Answer 1

[Mathieu NOE]

hello

most obviousy you didn't peak the right data in this plot (a copy paste from above that I fixed below, look for the comment : % <= here)

also i am not sure you need this . seems to me the output is better without

also , the parameters used with findpeaks seems not optimal as you don't pick the 4 dominant peaks .

    % Adjust the peak locations to account for the starting index of the response signal
    %Xpk = Xpk + 1;

clear

% Define the mass, damping, and stiffness matrices
load C.mat
load K.mat
load M.mat
% Define the input force as a sinusoidal function of frequency w in Hz
f = linspace(0, 40, 100); % linearly spaced frequency vector in Hz
w = 2*pi*f; % angular frequency vector in rad/s
F = @(t) [sin(2*pi*f(1)*t); zeros(3, 1)];
% Define the initial conditions
x0 = zeros(4, 1);
v0 = zeros(4, 1);
y0 = [x0; v0];
ind = 100;
% Define the time span
tspan = [0 20];
% Preallocate arrays to store the magnitude and phase of the response
mag = zeros(size(w));
phase = zeros(size(w));
X = zeros(length(w), 4);
% Loop over the frequency vector and compute the response at each frequency
for i = 1:length(w)
    % Define the input force at the current frequency
    F = @(t) [sin(w(i)*t); zeros(3, 1)];
    
    % Solve the differential equation to obtain the steady-state response
    [~, y] = ode45(@(t, y) vibration_equation(t, y, M, C, K, F), tspan, y0);
    x = y(end, 1:4)';
    
    
    % Compute the magnitude, phase, and displacement of the response at the current frequency
    mag(i) = abs(x(1));
    phase(i) = angle(x(1)) * 180/pi;
    X(i,:) = x';
    mpp = max(mag(1,:))/30;                                          
        [Ypk,Xpk,Wpk,Ppk] = findpeaks(mag(1,:),'MinPeakProminence',mpp, 'WidthReference','halfheight');  
        Xpk = Xpk + min(size(mag)) - 1;
        [~,mxidx] = maxk(Ppk,4); % to get top 4 peak prominences
        Ypkc{i} = Ypk(mxidx);     % peak amp                                                     
        Xpkc{i} = Xpk(mxidx);     % peak locations                                                       
        Wpkc{i} = Wpk(mxidx);     % peak widths                                                    
        Ppkc{i} = Ppk(mxidx);     % peak prominences                                                                                                  
        freq{i} = f(Xpkc{i});     % Frequency 
end
figure;
subplot(2, 1, 1);
%plot(f, mag);
plot(f,mag,f(Xpk),Ypk,'dr');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
ylim([0 3.5/10000]);
title('Frequency Response Function');
grid on;
subplot(2, 1, 2);
plot(f, mod(phase+180,360)-180,f(Xpk),Ypk,'dr');
xlabel('Frequency (Hz)');
ylabel('Phase (deg)');
ylim([-180, 180]);
grid on;
Ypkc = cell(1,4);
Xpkc = cell(1,4);
Wpkc = cell(1,4);
Ppkc = cell(1,4);
for i = 1:4
    % Find the minimum peak prominence as a fraction of the maximum response
    mpp(1,i) = max(X(:,i))/40;
    
    % Find the peaks in the response signal
    [Ypk,Xpk,Wpk,Ppk] = findpeaks(X(:,i), 'MinPeakProminence', mpp(1,i), 'WidthReference', 'halfheight');
    
    % Adjust the peak locations to account for the starting index of the response signal
    %Xpk = Xpk + 1;
    
    % Find the top 4 peak prominences and store the peak information
    [~, mxidx] = maxk(Ppk, 4);
    Ypkc{i} = Ypk(mxidx);
    Xpkc{i} = Xpk(mxidx);
    Wpkc{i} = Wpk(mxidx);
    Ppkc{i} = Ppk(mxidx);
end
% Plot the magnitude, phase, and displacement of the FRF in Hz
figure;
subplot(2, 2, 1);
plot(f, abs(X(:,1)), 'b');
hold on;
plot(f(Xpkc{1}), Ypkc{1}, 'dr'); % <= here
xlabel('Frequency (Hz)');
ylabel('Displacement (m)');
ylim([0 3.5/10000]);
title('X1');
grid on;
subplot(2, 2, 2);
plot(f, abs(X(:,2)),'b');
hold on;
plot(f(Xpkc{2}), Ypkc{2}, 'dr'); % <= here
xlabel('Frequency (Hz)');
ylabel('Displacement (m)');
ylim([0 3.5/10000]);
title('X2');
grid on;
subplot(2, 2, 3);
plot(f, abs(X(:,3)),'b');
hold on;
plot(f(Xpkc{3}), Ypkc{3},'dr'); % <= here
xlabel('Frequency (Hz)');
ylabel('Displacement (m)');
ylim([0 3.5/10000]);
title('X3');
grid on;
subplot(2, 2, 4);
plot(f, abs(X(:,4)),'b');
hold on;
plot(f(Xpkc{4}), Ypkc{4}, 'dr'); % <= here
xlabel('Frequency (Hz)');
ylabel('Displacement (m)');
ylim([0 3.5/10000]);
title('X4');
grid on;
clear C damp_ratio F i ind K K1_Th K2_Th M Natural_freq tspan mxidx

function dydt = vibration_equation(t, y, M, C, K, F)
% Compute the acceleration
x = y(1:4);
v = y(5:8);
a = M \ (F(t) - C*v - K*x);
% Compute the derivative of the state vector
dydt = [v; a];
end

[AL]

Thank you so much @Mathieu NOE for highlighting error and I will work on dominant peaks. Have a wonderful week, Once again I really appriciate your efforts.

[Mathieu NOE]

Glad I could be of some help

it was a minor problem

have a great day

[Mathieu NOE]

hello @AL

I could make some further improvements in the code

as you want explicitely the 4 dominant peaks , you can ask directly findpeaks to do so

[Ypk,Xpk,Wpk,Ppk] = findpeaks(mag(1,:),'NPeaks',4,'SortStr','descend');  

see below :

% Define the mass, damping, and stiffness matrices
load C.mat
load K.mat
load M.mat

% Define the input force as a sinusoidal function of frequency w in Hz
f = linspace(0, 40, 100); % linearly spaced frequency vector in Hz
w = 2*pi*f; % angular frequency vector in rad/s
F = @(t) [sin(2*pi*f(1)*t); zeros(3, 1)];

% Define the initial conditions
x0 = zeros(4, 1);
v0 = zeros(4, 1);
y0 = [x0; v0];
ind = 100;

% Define the time span
tspan = [0 20];

% Preallocate arrays to store the magnitude and phase of the response
mag = zeros(size(w));
phase = zeros(size(w));
X = zeros(length(w), 4);

% Loop over the frequency vector and compute the response at each frequency
for i = 1:length(w)
    % Define the input force at the current frequency
    F = @(t) [sin(w(i)*t); zeros(3, 1)];
    
    % Solve the differential equation to obtain the steady-state response
    [~, y] = ode45(@(t, y) vibration_equation(t, y, M, C, K, F), tspan, y0);
    x = y(end, 1:4)';
    
    % Compute the magnitude, phase, and displacement of the response at the current frequency
    mag(i) = abs(x(1));
    phase(i) = angle(x(1)) * 180/pi;
    X(i,:) = x';
    mpp = max(mag(1,:))/30;                                          

end

% find the 4 dominant peaks
[Ypk,Xpk,Wpk,Ppk] = findpeaks(mag(1,:),'NPeaks',4,'SortStr','descend');  
Xpk = Xpk + min(size(mag)) - 1;

figure;
subplot(2, 1, 1);
%plot(f, mag);
plot(f,mag,f(Xpk),Ypk,'dr');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
ylim([0 3.5/10000]);
title('Frequency Response Function');
grid on;
subplot(2, 1, 2);
plot(f, mod(phase+180,360)-180,f(Xpk),Ypk,'dr');
xlabel('Frequency (Hz)');
ylabel('Phase (deg)');
ylim([-180, 180]);
grid on;

% second plot with all four FRF's
Ypkc = cell(1,4);
Xpkc = cell(1,4);
Wpkc = cell(1,4);
Ppkc = cell(1,4);

for i = 1:4
    % Find the peaks in the response signal
    [Ypk,Xpk,Wpk,Ppk] = findpeaks(abs(X(:,i)),'NPeaks',4,'SortStr','descend');  
    Ypkc{i} = Ypk;
    Xpkc{i} = Xpk;
    Wpkc{i} = Wpk;
    Ppkc{i} = Ppk;
end

% Plot the magnitude, phase, and displacement of the FRF in Hz
figure;
subplot(2, 2, 1);
plot(f, abs(X(:,1)), 'b');
hold on;
plot(f(Xpkc{1}), Ypkc{1}, 'dr'); % <= here
xlabel('Frequency (Hz)');
ylabel('Displacement (m)');
ylim([0 3.5/10000]);
title('X1');
grid on;
subplot(2, 2, 2);
plot(f, abs(X(:,2)),'b');
hold on;
plot(f(Xpkc{2}), Ypkc{2}, 'dr'); % <= here
xlabel('Frequency (Hz)');
ylabel('Displacement (m)');
ylim([0 3.5/10000]);
title('X2');
grid on;
subplot(2, 2, 3);
plot(f, abs(X(:,3)),'b');
hold on;
plot(f(Xpkc{3}), Ypkc{3},'dr'); % <= here
xlabel('Frequency (Hz)');
ylabel('Displacement (m)');
ylim([0 3.5/10000]);
title('X3');
grid on;
subplot(2, 2, 4);
plot(f, abs(X(:,4)),'b');
hold on;
plot(f(Xpkc{4}), Ypkc{4}, 'dr'); % <= here
xlabel('Frequency (Hz)');
ylabel('Displacement (m)');
ylim([0 3.5/10000]);
title('X4');
grid on;
clear C damp_ratio F i ind K K1_Th K2_Th M Natural_freq tspan mxidx
function dydt = vibration_equation(t, y, M, C, K, F)
% Compute the acceleration
x = y(1:4);
v = y(5:8);
a = M \ (F(t) - C*v - K*x);
% Compute the derivative of the state vector
dydt = [v; a];
end