使用最小二乘法可以解決的問(wèn)題之
將一個(gè)方波分解為
asin(x)+bsin(2x)+...+csin(nx)+dcos(x)+ecos(2x)+...+fcos(nx)
產(chǎn)生一個(gè)方波：
使用周期函數(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，方便進(jìn)行矩陣的運(yùn)算
方波
將一個(gè)方波分解為
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
對(duì)比選擇不同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
可以看出n值越大，分解后的函數(shù)越接近方波函數(shù)
完整程序：
http://pan.baidu.com/s/1ckHTYu
提取密碼：kwrv
