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