Оптимизация окна Фурье за один цикл. Мой код неэффективен

Question

Оптимизация окна Фурье за один цикл. Мой код неэффективен

Добрый день

РЕДАКТИРОВАТЬ:Что я хочу: из любой формы волны тока / напряжения в системе питания (PS) мне нужны отфильтрованные (фундаментальные) значения RMS значений 50 Гц (и фактически их углы). Измеренный ток содержит все гармоники от 100 Гц до 1250 Гц в зависимости от оборудования. Невозможно анализировать и рассчитывать с помощью волны с этими гармониками, ваша ошибка становится настолько большой (в зависимости от оборудования), что защитное оборудование PS вычисляет неверные величины. Прилагаемый сигнал также содержит МНОГО многих других частотных компонентов.

Моя цель: реле защиты PS являются особенными и рассчитывают окно в 20 мс за очень короткое время. Я не пытаюсь этого понять. Я использую внешнюю технологию записи и проверяю, что реле видят правду и работают правильно. Таким образом, мне нужно делать то, что они делают, и сохранять только значения 50 Гц без каких-либо гармоник и постоянного тока.

Важный ожидаемый результат: учитывая любую частотную составляющую, которая МОЖЕТ быть в сигнале, я хочу увидеть величину любой данной гармоники (150,250 - амплитуды 3-й гармоники и 5-я гармоника основной гармоники), а также величину постоянного тока. Это скажет мне, какой тип оборудования PS может вводить эти частоты. Важно, чтобы я указывал частоту, и ответ был вектором этой частоты только со всеми другими значениями, отфильтрованными OUT. Среднеквадратичное значение основной гармоники и среднеквадратичное значение отличаются в 4000 А (только 50 Гц) и 4500 А (с другими включенными частотами)

Этот код вычисляет значение Фурье за один цикл (RMS) для заданной частоты. Я думаю, иногда его называют фильтром Фурье? Я использую его для анализа аналогов Power System 50Hz/0Hz/150Hz. (Ответы были протестированы и представляют собой правильные основные значения RMS. ( https://users.wpi.edu/~goulet/Matlab/overlap/trigfs.html)

Для большой выборки код очень медленный. Для 55000 точек данных требуется 12 секунд. Для 3 напряжений и 3 тока это будет ОЧЕНЬ медленно. Я просматриваю сотни записей в день.

Как мне его улучшить? Какие советы и хитрости / библиотеки Python есть для добавления моих списков / массива. (Также не стесняйтесь переписывать или использовать код). Я использую код для выделения определенных значений из сигнала в заданное время. (что похоже на чтение значений из специализированной программы для анализа энергосистемы)Отредактировано: с учетом того, как я загружаю файлы и использую их, код работает с их вставкой:

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()

Вставленный код работает без проблем, учитывая прикрепленные файлы. С уважением

РЕДАКТИРОВАТЬ: Здесь вы найдете тестовый файл csvfile и COMTRADE TEST:CSV: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

0

python signal-processing simpowersystems

Источник

user12307809 20 ноя '19 в 14:42

1 ответ

Решение

Другие вопросы по тегам python signal-processing simpowersystems

user3067485 20 ноя '19 в 19:08 2019-11-20 19:08 · Accepted Answer · 2019-11-20 19:08

Предисловие

Как я сказал в своем предыдущем комментарии:

Ваш код в основном полагается на forцикл с большим количеством операций индексации и скалярных операций. Вы уже импортировалиnumpy поэтому вам следует воспользоваться векторизацией.

Этот ответ - начало вашего решения.

Легкий 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

Или, альтернативно, мы можем вычислить его с помощью ДПФ (реализовано как rfft в numpy.fft):

# 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

Демонстрацию того, почему эта последняя формула работает, можно найти здесь. Это верно, поскольку из теоремы Парсеваля следует, что преобразование Фурье действительно сохраняет энергию.

Обе версии используют векторизованные функции, нет необходимости разделять действительную и мнимую части для выполнения вычислений, а затем повторной сборки в комплексное число.

MCVE: Окно

Я подозреваю, что вы хотите применить эту функцию в качестве окна для долгосрочного временного ряда, где значение 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'])

Скользящее среднеквадратичное значение вашего сигнала:

df['rms'] = df.rolling(int(fs/f)).agg(rms)

Периодически выбираемое RMS:

df['signal'].resample('1s').agg(rms)

Позднее возвращается:

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 Гц), простой решение может быть линейной detrend (с постоянной удалить и линейного тренда) с последующим Баттерворта фильтра (полосовой фильтр).

Мы генерируем синтетический сигнал с другими частотами и линейным трендом:

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.

Оптимизация окна Фурье за ​​один цикл. Мой код неэффективен

1 ответ

Предисловие

Легкий MCVE

MCVE: Окно

Преобразование сигнала

Оптимизация окна Фурье за один цикл. Мой код неэффективен