问傅里叶拟合系数

EN

import math
import numpy as np
#period of the signal
per=1.0
w = 2.0*np.pi/per
#number of fourier components.
nf = 5
fp = open("file.cat","r")
# m1 is the number of unknown coefficients.
m1 = 2*nf + 1
# Create empty matrices.
x = np.zeros((m1,m1))
y = np.zeros((m1,1))
xi = [0.0]*m1

# Read (time, value) from each line of the file.
for line in fp:
    t = float(line.split()[0])
    yi = float(line.split()[1])
    xi[0] = 1.0
    for k in range(1,nf+1):
        xi[2*k-1] = np.sin(k*w*t)
        xi[2*k] = np.cos(k*w*t)
    for j in range(m1):
        for k in range(m1):
            x[j,k] += xi[j]*xi[k]
        y[j] += yi*xi[j]
fp.close()
# Copy to big matrices.
X = np.mat( x.copy() )
Y = np.mat( y.copy() )
# Invert X and multiply by Y to get coefficients.
A = X.I*Y
A0 = A[0]
# Solution is A0 + Sum[ Amp*sin(k*wt + phi) ]
print "a[0] = %f" % A[0]
for k in range(1,nf+1):
    amp = math.sqrt(A[2*k-1]**2 + A[2*k]**2)
    phs = math.atan2(A[2*k],A[2*k-1])
    print "amp[%d] = %f phi = %f" % (k, amp, phs)

amp = np.array([np.sqrt((A[2*k-1])**2 + (A[2*k])**2) for k in range(1,nf+1)])
phs = np.array([math.atan2(A[2*k],A[2*k-1]) for k in range(1,nf+1)])

def term(t): return np.array([amp[k]*np.sin(k*w*t + phs[k]) for k in range(len(amp))])
Signal = np.array([A0+sum(term(phase[i])) for i in range(len(mag))])

def term(t): return np.array([amp[k]*np.sin((k+1)*w*t + phs[k]) for k in range(len(amp))])
plt.plot(phase,Signal,'r-',lw=3)