In [3]:

    

import numpy as np
from scipy.fftpack import fft
import matplotlib.pyplot as plt
%matplotlib


    

    




Using matplotlib backend: Qt5Agg

In [4]:

    

x = np.linspace(0, 1, 1400)
y = 7 * np.sin(2 * np.pi * 180 * x) + 2.8 *  np.sin(2 * np.pi * 390 * x) + 5.1 *  np.sin(2 * np.pi * 600 * x)


    

In [5]:

    

yy = fft(y)
yy


    

    
Out[5]:





array([  3.79962728e-12+0.j        ,   9.99966152e-05-0.0445618j ,
         4.00018214e-04-0.08913023j, ...,   9.00160110e-04+0.13371191j,
         4.00018214e-04+0.08913023j,   9.99966152e-05+0.0445618j ])

In [6]:

    

yy.real


    

    
Out[6]:





array([  3.79962728e-12,   9.99966152e-05,   4.00018214e-04, ...,
         9.00160110e-04,   4.00018214e-04,   9.99966152e-05])

In [8]:

    

yf = abs(yy)
yf


    

    
Out[8]:





array([  3.79962728e-12,   4.45619129e-02,   8.91311262e-02, ...,
         1.33714945e-01,   8.91311262e-02,   4.45619129e-02])

In [9]:

    

abs(1 + 1j)


    

    
Out[9]:





1.4142135623730951

In [10]:

    

plt.subplot(221)
plt.plot(x[0: 180], y[0:180])
plt.title('original wave')


    

    
Out[10]:





<matplotlib.text.Text at 0x1c34bb9b860>

In [11]:

    

yf1 = abs(yy) / len(x)
yf2 = yf1[range(int(len(x) / 2))]

xf = np.arange(len(x))
xf2=xf[range(int(len(x) / 2))]


    

In [12]:

    

plt.subplot(222)
plt.plot(xf,yf,'r')
plt.title('FFT of Mixed wave(two sides frequency range)',fontsize=7,color='#7A378B')  #注意这里的颜色可以查询颜色代码表


plt.subplot(223)
plt.plot(xf,yf1,'g')
plt.title('FFT of Mixed wave(normalization)',fontsize=9,color='r')


plt.subplot(224)
plt.plot(xf2,yf2,'b')
plt.title('FFT of Mixed wave)',fontsize=10,color='#F08080')


    

    
Out[12]:





<matplotlib.text.Text at 0x1c34c4ca2e8>

In [14]:

    

# http://blog.csdn.net/ouening/article/details/70339341  
# python实现傅立叶级数展开
x = np.mgrid[-10: 10.02: 0.02]
def fourier1():
    s = np.pi / 2
    print(s)

    for i in range(1,100,1):
        s0 = 2 / np.pi * (1 - (-1) ** i) / i ** 2 * np.cos(i * x) 
        s = s + s0
        if len(s) == 2:
            print('hehe')
            print(s)
#     plt.plot(x[0:10], s[0:10], 'orange', linewidth=0.6)
#     plt.title('fourier1')
#     plt.show()
#     print(s[0:10])
fourier1()


    

    




1.5707963267948966

In [ ]:

    




    

In [48]:

    

from pylab import *
from scipy.io import wavfile
import math

sampFreq, snd = wavfile.read('440_sine.wav')
print(sampFreq)
snd = snd / (2.**15)
# snd = snd

# snd.shape
s1 = snd[:, 0]

timeArray = arange(0, 5292.0, 1)   #[0s, 1s], 5292个点
timeArray = timeArray / sampFreq   #[0s, 0.114s]
timeArray = timeArray * 1000       #[0ms, 114ms]

plt.subplot(3, 1, 1)
plt.plot(timeArray, s1, color='k')
ylabel('Amplitude')
xlabel('Time (ms)')

n = len(s1)
p = fft(s1)         #执行傅立叶变换

nUniquePts = int(math.ceil((n + 1) / 2.0))
p = p[0: nUniquePts]
p = abs(p)

p = p / float(n)    #除以采样点数，去除幅度对信号长度或采样频率的依赖
p = p**2            #求平方得到能量

#乘2（详见技术手册）
#奇nfft排除奈奎斯特点
if n % 2 > 0:       #fft点数为奇
    p[1:len(p)] = p[1:len(p)]*2
else:               #fft点数为偶
    p[1:len(p)-1] = p[1:len(p)-1] * 2

freqArray = arange(0, nUniquePts, 1.0) * (sampFreq / n)
plt.subplot(3, 1, 2)
plt.plot(freqArray/1000, 10*log10(p), color='k')
xlabel('Freqency (kHz)')
ylabel('Power (dB)')

plt.subplot(3, 1, 3)
plt.plot(freqArray/1000, p, color='k')


    

    




44100






    
Out[48]:





[<matplotlib.lines.Line2D at 0x1c35d7a0400>]

In [ ]:

    




    

In [ ]: