Amplitude vs frequency curve
16 views (last 30 days)
Show older comments
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 Comments
Scott MacKenzie
on 2 Jul 2021
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);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/672238/image.jpeg)
Answers (0)
See Also
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!