单循环傅立叶 Window 优化。我的代码效率低下
One Cycle Fourier Window optimization. My code is inefficient
美好的一天
编辑：
我想要的：来自电力系统上的任何 current/voltage 波形（PS）我想要过滤后的 50Hz（基本）RMS 值幅度（以及有效的角度） .测量的电流包含从 100Hz 到 1250Hz 的所有谐波，具体取决于设备。无法使用具有这些谐波的波进行分析和计算，您的误差会变得如此之大（取决于设备），以至于 PS 保护设备计算出的数量不正确。附加的信号还涉及许多其他频率分量。
我的目标：PS 保护继电器很特殊，在很短的时间内计算出 20ms window。 I.m 不想得到这个。我正在使用外部记录技术并测试继电器看到的是真实的并且它们运行正常。因此，我需要做他们做的事，只保留 50Hz 值，没有任何谐波和直流。
重要的预期结果：给定信号中可能存在的任何频率分量，我想查看任何给定谐波的幅值（150,250 - 三次谐波幅值和五次谐波幅值）基本）以及DC的幅度。这将告诉我什么类型的 PS 设备可能注入这些频率。重要的是我提供了一个频率，答案是该频率的向量，所有其他值都被过滤掉了。
基波 RMS 与 RMS 的区别在于 4000A（仅 50Hz）和 4500A（包括其他频率）
此代码计算给定频率的单周期傅立叶值 (RMS)。我认为有时称为傅里叶滤波器？我将它用于 Power System 50Hz/0Hz/150Hz 模拟分析。 （答案已经过测试并且是正确的基本 RMS 值。(https://users.wpi.edu/~goulet/Matlab/overlap/trigfs.html)
对于大样本，代码非常慢。对于 55000 个数据点，需要 12 秒。对于 3 个电压和 3 个电流，这变得非常慢。我一天看100条记录。
如何增强它？有哪些 Python 提示和技巧/库可以附加我的 lists/array。
（也可以随意重写或使用代码）。我使用代码在给定时间从信号中挑选出某些值。 （这就像从电力系统分析的专门程序中读取值一样）
已编辑：通过我加载和使用文件的方式，代码可以粘贴它：
import matplotlib.pyplot as plt
import csv
import math
import numpy as np
import cmath
# FILES ATTACHED TO POST
filenamecfg = r"E:/Python_Practise/2019-10-21 13-54-38-482.CFG"
filename = r"E:/Python_Practise/2019-10-21 13-54-38-482.DAT"
t = []
IR = []
newIR=[]
with open(filenamecfg,'r') as csvfile1:
cfgfile = [row for row in csv.reader(csvfile1, delimiter=',')]
numberofchannels=int(np.array(cfgfile)[1][0])
scaleval = float(np.array(cfgfile)[3][5])
scalevalI = float(np.array(cfgfile)[8][5])
samplingfreq = float(np.array(cfgfile)[numberofchannels+4][0])
numsamples = int(np.array(cfgfile)[numberofchannels+4][1])
freq = float(np.array(cfgfile)[numberofchannels+2][0])
intsample = int(samplingfreq/freq)
#TODO neeeed to get number of samples and frequency and detect
#automatically
#scaleval = np.array(cfgfile)[3]
print('multiplier:',scaleval)
print('SampFrq:',samplingfreq)
print('NumSamples:',numsamples)
print('Freq:',freq)
with open(filename,'r') as csvfile:
plots = csv.reader(csvfile, delimiter=',')
for row in plots:
t.append(float(row[1])/1000000) #get time from us to s
IR.append(float(row[6]))
newIR = np.array(IR) * scalevalI
t = np.array(t)
def mag_and_theta_for_given_freq(f,IVsignal,Tsignal,samples): #samples are the sample window size you want to caclulate for (256 in my case)
# f in hertz, IVsignal, Tsignal in numpy.array
timegap = Tsignal[2]-Tsignal[3]
pi = math.pi
w = 2*pi*f
Xr = []
Xc = []
Cplx = []
mag = []
theta = []
#print("Calculating for frequency:",f)
for i in range(len(IVsignal)-samples):
newspan = range(i,i+samples)
timewindow = Tsignal[newspan]
#print("this is my time: ",timewindow)
Sig20ms = IVsignal[newspan]
N = len(Sig20ms) #get number of samples of my current Freq
RealI = np.multiply(Sig20ms, np.cos(w*timewindow)) #Get Real and Imaginary part of any signal for given frequency
ImagI = -1*np.multiply(Sig20ms, np.sin(w*timewindow)) #Filters and calculates 1 WINDOW RMS (root mean square value).
#calculate for whole signal and create a new vector. This is the RMS vector (used everywhere in power system analysis)
Xr.append((math.sqrt(2)/N)*sum(RealI)) ### TAKES SO MUCH TIME
Xc.append((math.sqrt(2)/N)*sum(ImagI)) ## these steps make RMS
Cplx.append(complex(Xr[i],Xc[i]))
mag.append(abs(Cplx[i]))
theta.append(np.angle(Cplx[i]))#th*180/pi # this can be used to get Degrees if necessary
#also for freq 0 (DC) id the offset is negative how do I return a negative to indicate this when i'm using MAGnitude or Absolute value
return Cplx,mag,theta #mag[:,1]#,theta # BUT THE MAGNITUDE WILL NEVER BE zero
myZ,magn,th = mag_and_theta_for_given_freq(freq,newIR,t,intsample)
plt.plot(newIR[0:30000],'b',linewidth=0.4)#, label='CFG has been loaded!')
plt.plot(magn[0:30000],'r',linewidth=1)
plt.show()
考虑到附加的文件，粘贴的代码运行顺利
问候
编辑：请在此处找到测试 csv 文件和 COMTRADE 测试文件：
格式文件：
https://drive.google.com/open?id=18zc4Ms_MtYAeTBm7tNQTcQkTnFWQ4LUu
商品贸易
https://drive.google.com/file/d/1j3mcBrljgerqIeJo7eiwWo9eDu_ocv9x/view?usp=sharing
https://drive.google.com/file/d/1pwYm2yj2x8sKYQUcw3dPy_a9GrqAgFtD/view?usp=sharing
前言
正如我在之前的评论中所说：
Your code mainly relies on a for loop with a lot of indexation and
scalar operations. You already have imported numpy so you should take
advantage of vectorization.
这个答案是您解决问题的开始。
轻量级 MCVE
首先我们为 MCVE 创建一个试验信号：
import numpy as np
# Synthetic signal sampler: 5s sampled as 400 Hz
fs = 400 # Hz
t = 5 # s
t = np.linspace(0, t, fs*t+1)
# Synthetic Signal: Amplitude is about 325V @50Hz
A = 325 # V
f = 50 # Hz
y = A*np.sin(2*f*np.pi*t) # V
然后我们可以使用通常的公式计算此信号的RMS：
# Actual definition of RMS:
yrms = np.sqrt(np.mean(y**2)) # 229.75 V
或者我们可以使用 numpy.fft 中的 DFT (implemented as rfft 来计算它：
# RMS using FFT:
Y = np.fft.rfft(y)/y.size
Yrms = np.sqrt(np.real(Y[0]**2 + np.sum(Y[1:]*np.conj(Y[1:]))/2)) # 229.64 V
可以找到为什么最后一个公式有效的演示here. This is valid because of the Parseval Theorem意味着傅里叶变换确实守恒。
两个版本都使用向量化函数，不需要将实部和虚部分开来进行计算，然后重新组合成一个复数。
MCVE：开窗
我怀疑您想将此函数作为 window 应用于 RMS 值即将发生变化的长期时间序列。然后我们可以使用提供时间序列商品的 pandas 库来解决这个问题。
import pandas as pd
我们封装了RMS函数：
def rms(y):
Y = 2*np.fft.rfft(y)/y.size
return np.sqrt(np.real(Y[0]**2 + np.sum(Y[1:]*np.conj(Y[1:]))/2))
我们生成阻尼信号（可变 RMS）
y = np.exp(-0.1*t)*A*np.sin(2*f*np.pi*t)
我们将试验信号包装到数据帧中以利用 rolling 或 resample 方法：
df = pd.DataFrame(y, index=t*pd.Timedelta('1s'), columns=['signal'])
您的信号的滚动 RMS 是：
df['rms'] = df.rolling(int(fs/f)).agg(rms)
周期性采样 RMS 为：
df['signal'].resample('1s').agg(rms)
后来的returns:
00:00:00 2.187840e+02
00:00:01 1.979639e+02
00:00:02 1.791252e+02
00:00:03 1.620792e+02
00:00:04 1.466553e+02
信号调节
为了满足仅保留基波谐波 (50 Hz) 的需求，一个简单的解决方案可能是线性 detrend (to remove constant and linear trend) followed by a Butterworth filter（带通滤波器）。
我们生成具有其他频率和线性趋势的合成信号：
y = np.exp(-0.1*t)*A*(np.sin(2*f*np.pi*t) \
+ 0.2*np.sin(8*f*np.pi*t) + 0.1*np.sin(16*f*np.pi*t)) \
+ A/20*t + A/100
然后我们调节信号：
from scipy import signal
yd = signal.detrend(y, type='linear')
filt = signal.butter(5, [40,60], btype='band', fs=fs, output='sos', analog=False)
yfilt = signal.sosfilt(filt, yd)
图形上它导致：
它恢复在 RMS 计算之前应用信号调节。
美好的一天
编辑： 我想要的：来自电力系统上的任何 current/voltage 波形（PS）我想要过滤后的 50Hz（基本）RMS 值幅度（以及有效的角度） .测量的电流包含从 100Hz 到 1250Hz 的所有谐波，具体取决于设备。无法使用具有这些谐波的波进行分析和计算，您的误差会变得如此之大（取决于设备），以至于 PS 保护设备计算出的数量不正确。附加的信号还涉及许多其他频率分量。
我的目标：PS 保护继电器很特殊，在很短的时间内计算出 20ms window。 I.m 不想得到这个。我正在使用外部记录技术并测试继电器看到的是真实的并且它们运行正常。因此，我需要做他们做的事，只保留 50Hz 值，没有任何谐波和直流。
重要的预期结果：给定信号中可能存在的任何频率分量，我想查看任何给定谐波的幅值（150,250 - 三次谐波幅值和五次谐波幅值）基本）以及DC的幅度。这将告诉我什么类型的 PS 设备可能注入这些频率。重要的是我提供了一个频率，答案是该频率的向量，所有其他值都被过滤掉了。 基波 RMS 与 RMS 的区别在于 4000A（仅 50Hz）和 4500A（包括其他频率）
此代码计算给定频率的单周期傅立叶值 (RMS)。我认为有时称为傅里叶滤波器？我将它用于 Power System 50Hz/0Hz/150Hz 模拟分析。 （答案已经过测试并且是正确的基本 RMS 值。(https://users.wpi.edu/~goulet/Matlab/overlap/trigfs.html)
对于大样本，代码非常慢。对于 55000 个数据点，需要 12 秒。对于 3 个电压和 3 个电流，这变得非常慢。我一天看100条记录。
如何增强它？有哪些 Python 提示和技巧/库可以附加我的 lists/array。 （也可以随意重写或使用代码）。我使用代码在给定时间从信号中挑选出某些值。 （这就像从电力系统分析的专门程序中读取值一样） 已编辑：通过我加载和使用文件的方式，代码可以粘贴它：
import matplotlib.pyplot as plt
import csv
import math
import numpy as np
import cmath
# FILES ATTACHED TO POST
filenamecfg = r"E:/Python_Practise/2019-10-21 13-54-38-482.CFG"
filename = r"E:/Python_Practise/2019-10-21 13-54-38-482.DAT"
t = []
IR = []
newIR=[]
with open(filenamecfg,'r') as csvfile1:
cfgfile = [row for row in csv.reader(csvfile1, delimiter=',')]
numberofchannels=int(np.array(cfgfile)[1][0])
scaleval = float(np.array(cfgfile)[3][5])
scalevalI = float(np.array(cfgfile)[8][5])
samplingfreq = float(np.array(cfgfile)[numberofchannels+4][0])
numsamples = int(np.array(cfgfile)[numberofchannels+4][1])
freq = float(np.array(cfgfile)[numberofchannels+2][0])
intsample = int(samplingfreq/freq)
#TODO neeeed to get number of samples and frequency and detect
#automatically
#scaleval = np.array(cfgfile)[3]
print('multiplier:',scaleval)
print('SampFrq:',samplingfreq)
print('NumSamples:',numsamples)
print('Freq:',freq)
with open(filename,'r') as csvfile:
plots = csv.reader(csvfile, delimiter=',')
for row in plots:
t.append(float(row[1])/1000000) #get time from us to s
IR.append(float(row[6]))
newIR = np.array(IR) * scalevalI
t = np.array(t)
def mag_and_theta_for_given_freq(f,IVsignal,Tsignal,samples): #samples are the sample window size you want to caclulate for (256 in my case)
# f in hertz, IVsignal, Tsignal in numpy.array
timegap = Tsignal[2]-Tsignal[3]
pi = math.pi
w = 2*pi*f
Xr = []
Xc = []
Cplx = []
mag = []
theta = []
#print("Calculating for frequency:",f)
for i in range(len(IVsignal)-samples):
newspan = range(i,i+samples)
timewindow = Tsignal[newspan]
#print("this is my time: ",timewindow)
Sig20ms = IVsignal[newspan]
N = len(Sig20ms) #get number of samples of my current Freq
RealI = np.multiply(Sig20ms, np.cos(w*timewindow)) #Get Real and Imaginary part of any signal for given frequency
ImagI = -1*np.multiply(Sig20ms, np.sin(w*timewindow)) #Filters and calculates 1 WINDOW RMS (root mean square value).
#calculate for whole signal and create a new vector. This is the RMS vector (used everywhere in power system analysis)
Xr.append((math.sqrt(2)/N)*sum(RealI)) ### TAKES SO MUCH TIME
Xc.append((math.sqrt(2)/N)*sum(ImagI)) ## these steps make RMS
Cplx.append(complex(Xr[i],Xc[i]))
mag.append(abs(Cplx[i]))
theta.append(np.angle(Cplx[i]))#th*180/pi # this can be used to get Degrees if necessary
#also for freq 0 (DC) id the offset is negative how do I return a negative to indicate this when i'm using MAGnitude or Absolute value
return Cplx,mag,theta #mag[:,1]#,theta # BUT THE MAGNITUDE WILL NEVER BE zero
myZ,magn,th = mag_and_theta_for_given_freq(freq,newIR,t,intsample)
plt.plot(newIR[0:30000],'b',linewidth=0.4)#, label='CFG has been loaded!')
plt.plot(magn[0:30000],'r',linewidth=1)
plt.show()
考虑到附加的文件，粘贴的代码运行顺利 问候
编辑：请在此处找到测试 csv 文件和 COMTRADE 测试文件： 格式文件： https://drive.google.com/open?id=18zc4Ms_MtYAeTBm7tNQTcQkTnFWQ4LUu
商品贸易 https://drive.google.com/file/d/1j3mcBrljgerqIeJo7eiwWo9eDu_ocv9x/view?usp=sharing https://drive.google.com/file/d/1pwYm2yj2x8sKYQUcw3dPy_a9GrqAgFtD/view?usp=sharing
前言
正如我在之前的评论中所说：
Your code mainly relies on a
forloop with a lot of indexation and scalar operations. You already have importednumpyso you should take advantage of vectorization.
这个答案是您解决问题的开始。
轻量级 MCVE
首先我们为 MCVE 创建一个试验信号：
import numpy as np
# Synthetic signal sampler: 5s sampled as 400 Hz
fs = 400 # Hz
t = 5 # s
t = np.linspace(0, t, fs*t+1)
# Synthetic Signal: Amplitude is about 325V @50Hz
A = 325 # V
f = 50 # Hz
y = A*np.sin(2*f*np.pi*t) # V
然后我们可以使用通常的公式计算此信号的RMS：
# Actual definition of RMS:
yrms = np.sqrt(np.mean(y**2)) # 229.75 V
或者我们可以使用 numpy.fft 中的 DFT (implemented as rfft 来计算它：
# RMS using FFT:
Y = np.fft.rfft(y)/y.size
Yrms = np.sqrt(np.real(Y[0]**2 + np.sum(Y[1:]*np.conj(Y[1:]))/2)) # 229.64 V
可以找到为什么最后一个公式有效的演示here. This is valid because of the Parseval Theorem意味着傅里叶变换确实守恒。
两个版本都使用向量化函数，不需要将实部和虚部分开来进行计算，然后重新组合成一个复数。
MCVE：开窗
我怀疑您想将此函数作为 window 应用于 RMS 值即将发生变化的长期时间序列。然后我们可以使用提供时间序列商品的 pandas 库来解决这个问题。
import pandas as pd
我们封装了RMS函数：
def rms(y):
Y = 2*np.fft.rfft(y)/y.size
return np.sqrt(np.real(Y[0]**2 + np.sum(Y[1:]*np.conj(Y[1:]))/2))
我们生成阻尼信号（可变 RMS）
y = np.exp(-0.1*t)*A*np.sin(2*f*np.pi*t)
我们将试验信号包装到数据帧中以利用 rolling 或 resample 方法：
df = pd.DataFrame(y, index=t*pd.Timedelta('1s'), columns=['signal'])
您的信号的滚动 RMS 是：
df['rms'] = df.rolling(int(fs/f)).agg(rms)
周期性采样 RMS 为：
df['signal'].resample('1s').agg(rms)
后来的returns:
00:00:00 2.187840e+02
00:00:01 1.979639e+02
00:00:02 1.791252e+02
00:00:03 1.620792e+02
00:00:04 1.466553e+02
信号调节
为了满足仅保留基波谐波 (50 Hz) 的需求，一个简单的解决方案可能是线性 detrend (to remove constant and linear trend) followed by a Butterworth filter（带通滤波器）。
我们生成具有其他频率和线性趋势的合成信号：
y = np.exp(-0.1*t)*A*(np.sin(2*f*np.pi*t) \
+ 0.2*np.sin(8*f*np.pi*t) + 0.1*np.sin(16*f*np.pi*t)) \
+ A/20*t + A/100
然后我们调节信号：
from scipy import signal
yd = signal.detrend(y, type='linear')
filt = signal.butter(5, [40,60], btype='band', fs=fs, output='sos', analog=False)
yfilt = signal.sosfilt(filt, yd)
图形上它导致：
它恢复在 RMS 计算之前应用信号调节。