50HZ陷波器python源码

#50HZ陷波器
#分子
numerator =[0.96897, -1.84310, 0.96897 ]

#分母
denominator =[1, -1.84310 ,0.93795]

import numpy as np
from matplotlib import pyplot as plt
import random

W_50_Hz = 2np.pi50 #定义频率W

#x_row_length = 2*np.pi/W_50_Hz; #2个周期 T=2π/w

one_cycle_sample_point=1000 #1个周期内采样1000个点
cycle = 5

total_sample_number=(one_cycle_sample_point*cycle)

x_row_length = cycle*0.02 #50HZ的周期是1/50=0.02

total_num_length = int(cycle*50000) #原始信号一共的点数

#1个周期内采样1000个点，采样间隔是
sample_interval=int(total_num_length/cycle/1000)
print("sample_interval is ",sample_interval)

#sample_fre = 1000#采样频率。

Sin50_data = [0,0]
#Sin70_data = [0,0]

noise = []
t = [0]
sample_point = [0,0]
sample_Sin50_data = [0,0]

filter_Sin50_data = [0,0]
#filter_Sin70_data = [0,0]

#W_70_Hz = 2np.pi70 #定义频率W

t=np.linspace(0,x_row_length,(total_num_length))#原始信号的时间间隔 x=linspace（x0，xn，m）表示在x0和xn之间等间隔取m个数

#参数loc(float)：正态分布的均值，对应着这个分布的中心。
#参数scale(float)：正态分布的标准差，对应分布的宽度，scale越大，正态分布的曲线越矮胖，scale越小，曲线越高瘦。
#参数size(int 或者整数元组)：输出的值赋在shape里，默认为None。
noise=np.random.normal(0,0.01,(total_num_length))#产生噪声

Sin50_data = np.sin(W_50_Hz*t) + noise#原始信号添加噪声

#sample_num=int(total_num_length/sample_fre)
#length = len(sample_Sin50_data)

index=0
for n in range(2,total_sample_number):
time_point=2np.pi50*n/one_cycle_sample_point
sample_Sin50_data.append( np.sin(time_point))
#sample_Sin50_data[index]=Sin50_data[n]
sample_point.append(time_point)
#sample_point[index]=t[n]
#index=index+1

sample_length=len(sample_point)
print(“sample point length is %d sample data length is”,len(sample_point),len(sample_Sin50_data))

for n in range(2,total_sample_number):# range（a,b,c） ,a为循环开始的数字(可不填，默认为0)，b为循环结束的后一位（c为正数时）的数字
now_data = 0.974025sample_Sin50_data[n] - 1.852705sample_Sin50_data[n-1] + 0.974025sample_Sin50_data[n-2] + 1.852705filter_Sin50_data[n-1] - 0.948050*filter_Sin50_data[n-2]
filter_Sin50_data.append(now_data)
#filter_Sin50_data[index]=now_data
#index=index+1

#del filter_Sin50_data[0 : 2]

print(“sample point length is\n\r”,len(sample_point))
print(“filter_Sin50_data length is\n\r”,len(filter_Sin50_data))
plt.plot(sample_point, sample_Sin50_data,color=‘red’, label=“Sin50_data”)
plt.plot(sample_point ,filter_Sin50_data,color=‘blue’, label=“filter_Sin50_data”)

plt.legend()

plt.show()