-
Notifications
You must be signed in to change notification settings - Fork 3
/
Frequency_spectrum.m
48 lines (43 loc) · 2.07 KB
/
Frequency_spectrum.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
function [] = Frequency_spectrum(a, b, c, force1, force2, irr_freq, ini_cond)
% Find frequency spectrum of solution generated by given parameters and initial condition.
%
% Inputs:
% a = -1*Damping strength
% b = Linear stiffness
% c = non-linearity in restoring force
% force1 = amplitude of force in rational frequency driving force (2*pi Hz)
% force2 = amplitude of force in irrational frequecy driving force (2*pi*irr_freq) Hz)
% irr_freq = value of irrational driving frequency
% init_cond = initial conditions for dynamical system
%
% Outputs:
% None.
%
% Figures:
% Plot of frequency spectrum of solution (amplitude against frequency)
%
% To display the frequency spectrum, first the data itself is generated by the Duffing_solution
% function with the paramters input here. Then we perform the discrete Fourier transform using
% Matlab's FFT algorithm on only the last 90% of the data (to remove transients), shift the data
% into a spectrum reflected about x = 0 and finally we double and divde the amplitudes by number of
% points. To get the x-axis right, the radians per second in the wave data given have to be found
% and taken into account. The shifted specturm is chosen because it literally just looks better, and
% the axes are restriced here to make the important low-frequency peaks more visible.
% Initialising values needed
[pos, ~] = Duffing_solution(a, b, c, force1, force2, irr_freq, ini_cond);
num_points = length(pos);
% Calculating the DFT, and fftshifting it
amplitudes = fft(pos(ceil(0.9*num_points) : num_points));
amplitudes = fftshift(amplitudes);
amplitudes = abs(2*amplitudes/num_points);
% Setting up the X-axis (frequency)
num_points = length(amplitudes);
rads_per_sec = 100;
frequencies = linspace(-num_points/2, num_points/2, num_points) * (rads_per_sec/num_points);
plot(frequencies, amplitudes);
xlabel('Frequency(Hz)');
ylabel('Amplitude');
title('Frequency spectrum of solution to duffing oscillator');
% Restrict view to make important part clearer
axis([-3, 3, 0, max(amplitudes)*1.2]);
end