用Python实现小波变换和傅里叶变换消除白噪声
算法原理
(1)小波变换
(2)傅里叶变换
DFT的正变换计算公式为:
DFT的逆变换计算公式为:
指数部分根据 Euler公式计算:
# -*- coding: utf-8 -*-
"""
Created on Wed Mar 18 16:16:28 2020
@author: 侯明会
"""
import math
import matplotlib.pyplot as plt
N=128
data=[0]*N
#在文件中读取128个数据
fo=open("lorenza.dat","r+")
list_x=list(fo.read().splitlines(False))
fo.close()
for i in range(N):
data[i]=float(list_x[i])
#小波变换
def DWT():
k=3 #k值设为3
A=[[]]*k
D=[[]]*k
print("数据图为:")
plt.plot(data)
plt.axis([0,128,0,25])
plt.show()
#小波分解
for i in range(int(N/2)):
A[0]=A[0]+[(data[2*i]+data[2*i+1])/math.sqrt(2)]
D[0]=D[0]+[(data[2*i]-data[2*i+1])/math.sqrt(2)]
n=N/2
for j in range(1,k):
for i in range(int(n/2)):
A[j]=A[j]+[(A[j-1][2*i]+A[j-1][2*i+1])/math.sqrt(2)]
D[j]=D[j]+[(A[j-1][2*i]-A[j-1][2*i+1])/math.sqrt(2)]
n=n/2
DWT_x=A[j]
for i in range(k):
DWT_x=DWT_x+D[i]
print("分解图为:")
plt.plot(DWT_x)
plt.axis([0,128,-4,60])
plt.show()
#滤波
Deltaw=10
for i in range(N):
if DWT_x[i]<Deltaw:#把DWT_x中小于Deltaw的数值归零
DWT_x[i]=0
print("滤波后图为(Deltaw=10):")
plt.plot(DWT_x)
plt.axis([0,128,-4,60])
plt.show()
#逆变换
for i in range(len(A[j])):#把A[j]中小于Deltaw的数值归零
if A[j][i]<Deltaw:
A[j][i]=0
for i in range(k):
for l in range(len(D[i])):#把D[k]中小于Deltaw的数值归零
if D[i][l]<Deltaw:
D[i][l]=0
for m in range(k-1,-1,-1):
for i in range(len(A[m])):
if m-1>-1:
A[m-1][2*i]=(A[m][i]+D[m][i])/math.sqrt(2)
A[m-1][2*i+1]=(A[m][i]-D[m][i])/math.sqrt(2)
for i in range(len(A[1])):
A[0][2*i]=(A[1][i]+D[1][i])/math.sqrt(2)
A[0][2*i+1]=(A[1][i]-D[1][i])/math.sqrt(2)
new_data=[0]*N
for i in range(len(A[0])):
new_data[2*i]=(A[0][i]+D[0][i])/math.sqrt(2)
new_data[2*i+1]=(A[0][i]-D[0][i])/math.sqrt(2)
print("逆变换后图为:")
plt.plot(new_data)
plt.axis([0,128,0,25])
plt.show()
return
#傅里叶变换
def DFT():
X=[0]*N
FF=[]
print("数据图为:")
plt.plot(data)
plt.axis([0,128,0,25])
plt.show()
for k in range(N):
F=[[],[]]#一个存实数 一个存虚数
for i in range(N):
a=-k*i*2*math.pi/N
F[0]=F[0]+[math.cos(a)*data[i]]
F[1]=F[1]+[math.sin(a)*data[i]]
f1=f2=0
for j in range(N):
f1=f1+F[0][j]
f2=f2+F[1][j]
FF.append([f1,f2])
X[k]=math.sqrt(f1**2+f2**2)
print("正变换后图为:")
plt.plot(X)
plt.show()
#滤波
Deltaf = 100
for k in range(N):
if X[k]<Deltaf:#把X中小于Deltaf的数值归零
X[k]=0
FF[k]=[0,0]
print("滤波后图为(Deltaf = 100):")
plt.plot(X)
plt.show()
#逆变换
X1=[0]*N
for n in range(N):
F1=[[],[]]#一个存实数 一个存虚数
for i in range(N):
a=n*i*2*math.pi/N
F1[0]=F1[0]+[(FF[i][0]*math.cos(a)-FF[i][1]*math.sin(a))/N]
F1[1]=F1[1]+[(FF[i][1]*math.cos(a)+FF[i][0]*math.sin(a))/N]
f11=f12=0
for j in range(N):
f11=f11+F1[0][j]
f12=f12+F1[1][j]
X1[n]=math.sqrt(f11**2+f12**2)
print("逆变换后图为:")
plt.plot(X1)
plt.axis([0,128,0,25])
plt.show()
return
print("小波变换过程图如下:")
DWT()
print("傅里叶变换过程图如下:")
DFT()