58 lines
1.3 KiB
Python
58 lines
1.3 KiB
Python
import numpy as np
|
|
|
|
|
|
def fft(x, fs, X_nom=None, epsilon=None):
|
|
|
|
# Implementation from IEC 61000-4-7:2002/AMD1:2008
|
|
|
|
|
|
# Data length
|
|
N = len(x)
|
|
# Number of positive-frequency bins
|
|
K = int(np.floor(N/2))
|
|
# Frequency axis (0 .. fs/2)
|
|
#freq = np.linspace(start=0, stop=fs/2, num=K+1, endpoint=True)
|
|
freq = np.arange(K + 1) * fs / N
|
|
|
|
# allocate
|
|
a = np.zeros(K + 1)
|
|
b = np.zeros(K + 1)
|
|
c = np.zeros(K + 1)
|
|
Y_C = np.zeros(K + 1)
|
|
phi = np.zeros(K + 1)
|
|
|
|
n = np.arange(N)
|
|
|
|
# DC
|
|
c[0] = np.mean(x) # c0 per IEC
|
|
a[0] = 2 * c[0] # not really used; just for completeness
|
|
b[0] = 0.0
|
|
|
|
# k = 1..K
|
|
for k in range(1, K+1):
|
|
angle = 2 * np.pi * k * n / N
|
|
a[k] = (2/N) * np.sum(x * np.cos(angle))
|
|
b[k] = (2/N) * np.sum(x * np.sin(angle))
|
|
|
|
# Nyquist (if N even): do NOT apply the 2/N doubling
|
|
if (N % 2 == 0) and (k == K):
|
|
a[k] *= 0.5
|
|
b[k] *= 0.5
|
|
c[k] = np.sqrt(a[k]*a[k] + b[k]*b[k])
|
|
|
|
# RMS value calculated in Eq.2.
|
|
Y_C[k] = c[k] / np.sqrt(2)
|
|
|
|
|
|
# Phase: apply dead-band if provided, else always compute
|
|
if (X_nom is not None) and (epsilon is not None):
|
|
if (np.abs(a[k]) <= epsilon * X_nom) and (np.abs(b[k]) <= epsilon * X_nom):
|
|
phi[k] = 0.0
|
|
continue
|
|
|
|
# IEC quadrant handling (equivalent to the piecewise definition)
|
|
phi[k] = np.arctan2(a[k], b[k])
|
|
|
|
|
|
|
|
return freq, a, b, c, Y_C, phi |