使用最小二乘法可以解決的問題之
將一個方波分解為
asin(x)+bsin(2x)+...+csin(nx)+dcos(x)+ecos(2x)+...+fcos(nx)

產(chǎn)生一個方波：
使用周期函數(shù)sin產(chǎn)生方波，sinx>0,y=-1,sinx<0,y=1

x = np.linspace(-10,10,300)
y=[]
for i in x:
    if np.sin(i)>0:#調(diào)用sin，cos要使用np.sin，np.cos
        y.append(-1)
    else:
        y.append(1)
y=np.array(y)#需要把list轉(zhuǎn)化成array，方便進行矩陣的運算

方波

將一個方波分解為
asin(nx),bcon(nx)的線性組合，
如：
asin(x)+bsin(2x)+...+csin(nx)+dcos(x)+ecos(2x)+...+fcos(nx)
要求的是系數(shù)a，b...c,d...e,f組成的矩陣，輸入量是方波（x，y）與n的值。
所以定義函數(shù)：

def fourier(x,y,n):
    return ym
#返回值為asin(nx),bcon(nx)的線性組合，
#即，系數(shù)a，b...c,d...e,f組成的矩陣與x1（sin(nx),con(nx)）的乘積

函數(shù)fourier

def fourier(x,y,n):
    x1=[]#（sin(nx),con(nx)）
    for i in xrange(n):
        x1.append(np.sin(x*i+x))
        x1.append(np.cos(x*i+x))
    m=np.mat(x1).T#使用np.mat方便矩陣的連乘
    y.shape=(y.shape[0],1)
    p=m*np.linalg.inv(m.T*m)*m.T*y
    ym=np.array(p)#將矩陣轉(zhuǎn)換成array，與前面統(tǒng)一
    ym.shape=(ym.shape[0],)
    return ym

對比選擇不同n值得分解結(jié)果：

plt.plot(x,y,color="g",label=u'方波')
plt.plot(x,fourier(x,y,3),color='r',label='3')
plt.plot(x,fourier(x,y,8),color='b',label='8')
plt.plot(x,fourier(x,y,23),color='k',label='23')
plt.legend()
plt.axis('equal')
plt.show()

Paste_Image.png