Amplitude vs frequency curve
13 次查看(过去 30 天)
显示 更早的评论
Hey guys ! I would like your help to plot this frequency response function.
Just so you guys understand, I'm trying to plot the frequency response curve through ode45. I have a one degree of freedom system with its constants and I'm driving this system with a sinusoidal force that will do a frequency sweep.
I am sending the analytic result, as for my code attempt. The curve using the fft appears to be coherent, but the y axis has a very different unit.
function simulation
clear all
close all
clc
m = 0.086; % kg
k = 166.3629; % N/m
c = 0.1664 ; % N.s/m
F = 0.6385;% N
% Time
Fs=400;
tspan =0:1/Fs:250-1/Fs;
f0 = 4; % Hz
f1 = 9; % Hz
% analytical solution
w = f0:0.2:f1; %Hz
d1 = (k - m*(w*2*pi).^2).^2 + (c*w*2*pi).^2;
X = F./sqrt(d1);
% IC
x0 = 0; v0 = 0;
IC2 = [x0;v0];
% numerical integration
[time2,state_values2] = ode45(@h,tspan,IC2);
x = state_values2(:,1);
figure(1)
plot(time2,x),xlabel('time(s)'),ylabel('displacement(m)')
figure(2)
n=ceil(log2(length(x)));
fx=fft(x,2^n);
fx=2*fx/length(x); % This operation is Adjusting the Magnitudes
f=(Fs/2^n)*(0:2^(n-1)-1);
plot(f,abs(fx(1:2^(n-1))),w,X,'-.k'), xlabel(' Frequency (Hz)'), ylabel(' X (m)');
l1 = ' using fft';
l2 = ' analytical solution';
legend(l1,l2);
xlim([f0 f1])
end
function sdot1 = h(t,x)
m = 0.086; % kg
k = 166.3629; % N/m
c = 0.1664 ; % N.s/m
f0 = 4; % Hz
f1 = 9; % Hz
F = 0.6385;% N
a = (f1 - f0)/250;
sdot1 = [x(2);
(F.*sin(2*pi*(a*t/2 + f0)*t) - c.*x(2) - k*x(1))/m];
end
4 个评论
Scott MacKenzie
2021-7-2
Is the y-axis magnitude really that important? Both solutions identify the signal at 7 Hz. Isn't that the key result? Sorry, if I'm completely off base here; this view is just based on my limited experience with spectrum analyses. I also notice that you resized the magnitudes after applying the fft. If you do so before, the result is closer to what you are looking for. Good luck.
x = rescale(x,-1, 1);
fx=fft(x,2^n);
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!