The notebook by defult uses the .png format for the inline plots, but if you've not that many points in your plots the SVG format can be prettier.
In [1]:
%config InlineBackend.figure_format = 'svg'
Our signal, or loading, is 3s long, but we need a stretch of zeros to damp out the response and simulate rest initial conditions, so we use a length, hence also a period T, of 8s.
The number of samples per second is here given as 512, note that the signal length in terms of samples is a power of 2 (here $2^{12}$).
In [2]:
T = 8
sps = 512
The array t contains the times at which our signal was sampled, the load is computed as a list by a procedure peculiar to computer languages Haskell and Python, known as list comprehension.
In [3]:
t = arange(0,T,1./sps)
p = [0 if x>3 else (60000.-20000.*x if x>1 else x*40000.) for x in t]
Am I sure that the list p contains the values of the loading?
Let's try to plot p vs t...
In [4]:
plot(t, p) ; xlabel("t/s") ; ylabel("p(t)/N")
print len(t), len(p)
It looks OK...
Now, the fast Fourier transform of the sequence p is computed, and given a name, P.
It is customary to denote Fourier pairs by the same letter, the small letter for the time domain representation and the capital letter for the frequency domain representation.
In [5]:
P = fft.fft(p)
iP = fft.ifft(P)
I have computed also the inverse FFT of the FFT of the loading, naming it iP, it is a sequence of complex numbers and here we plot the real and the imaginary part of each component versus time.
In [6]:
plot(t,real(iP),'k',t,imag(iP),'y') ; xlabel("t/s") ; ylabel("p(t)/N") ;
It seems OK...
Next, we use a convenience function to compute a sequence of frequencies (in Hertz!) associated with the components of P, the FFT of p. The parameters are the number of points and the sampling interval..
In [7]:
f = fft.fftfreq(T*sps, 1./sps)
The x axis is streching over the interval $-f_\text{Ny}$, $+f_\text{Ny}$
In [8]:
plot(f,real(P),'k', f,imag(P),'b')
xlim(-256,256) ; xticks(range(-256,257,64))
xlabel("f/Hz") ; ylabel("P(f)") ; grid() ;
The plot above is not much clear, in the next 3 plots we zoom near the origin of the frequency axis to have a bit more of detail. There are 3 plots, first the absolute value of P vs f, then the real part and finally the imaginary part of P, versus f.
In [9]:
plot(f,abs(P),'k')
xticks(range(-4,5,2))
xlabel("f/Hz") ; ylabel("abs(P(f))")
xlim(-4,4) ; ylim(0,3.3E7)
grid() ;
In [10]:
plot(f,real(P),'k')
xticks(range(-4,5,2))
xlabel("f/Hz") ; ylabel("real(P(f))")
xlim(-4,4) ; ylim(-3.3E7,3.3E7)
grid() ;
In [11]:
plot(f,imag(P),'k')
xticks(range(-4,5,2))
xlabel("f/Hz") ; ylabel("imag(P(f))")
xlim(-4,4) ; ylim(-3.3E7,3.3E7)
grid() ;
Until now, we did without the SDOF, now it's time to describe it and derive its response function.
All the parameters are the same as in the excel example, we compute k because we need it to normalize the response.
In [12]:
z = 0.1; fn = 1/0.6 ; m =6E5 ; wn = fn*2*pi ; k = m*wn**2
def H(f):
b = f/fn
return 1./((1-b*b)+1j*(2*z*b))
As usual, we plot the response function, or rather the absolute value of, against a short span of the frequency axis, centered about the origin, to show the details of the response function itself.
In [13]:
plot(f,abs(H(f))) ; xlabel("f/Hz") ; xlim(-8,8) ; ylabel("H(f)") ;
The FFT of the response is computed multiplying, term by term, P by the response or transfer function (again, that's done using a list comprehension), then we compute the IFFT of X to obtain x, the representation of the periodic response in the time domain.
In [14]:
X = [ P_*H(f_) for P_, f_ in zip(P,f)]
x = ifft(X)/k
In the end, we remain with the task of plotting the response function, that is the real part of x. Just to be certain we plot also the imaginary part of x, multiplied by 10**13, so we can be sure that it is negligible with respect to the real part
In [15]:
plot(t,real(x), t,imag(x)*1E13)
grid(); xlabel("t/s"); ylabel("x/m");
In [ ]: