The modification which you have suggested are working perfectly and I really appricate youe efforts. Have a wonderful day.
-AL.
finding peak points while signal processing
4 次查看（过去 30 天）
显示 更早的评论
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.
采纳的回答
Mathieu NOE
2023-3-20
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
3 个评论
Mathieu NOE
2023-3-20
移动：Mathieu NOE
2023-3-20
Glad I could be of some help
it was a minor problem
have a great day
Mathieu NOE
2023-3-20
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
更多回答（0 个）
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Pulsed Waveforms 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
您也可以从以下列表中选择网站：
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)