Frequency Domain Filtering Fundamentals
Introduction
Filtering in the frequency domain consists of modifying the Fourier transform of an signal (can be a image, a media file, a light curve…) and then taking the inverse tranform to obtained the filtered result. Thus, given a digital signal, \(f(x)\), of length \(M\), the basic filtering equation is:
where \(\mathcal{F}^{-1}\) is the Inverse discrete Fourier transform, \(F(u)\) is the Fourier transform of the given function \(f(x)\) (input), \(H(u)\) is the filter transfer function, and \(g(x)\) is the filtered signal (output). Both \(F\), \(H\) and \(G\) are arrays of size \(M\), the same as the input signal. The product \(H(u)F(u)\) is formed using array multiplication; that is, \(G(i, k)=H(i, k)F(i, k)\).
Expand LightCurves borders
To perform the filtering of the light curves in this work, we will use the Fast Fourier Transform (FFT) algorithm in order to convert the signal in its original domain in the frequency domain. The choice of FFT is simple, this algorithm has complexity of \(O(N log\ N)\), while the Discrete-time Fourier transform has complexity \(O(N^2)\), where \(N\) is the input data size, it means, FFT is much more computationally efficient.
However, the FFT has some peculiarities that can lead to problemas and unexpected results. The algorithm is based on the sucessive folding method [Gonzales and Woods, 1992], in other words, it assumes that the input signal is periodic and that the number of samples is an integer and a power of 2, but the second premise is not a requirement general. The first conjecture is quite relevant to not having distortion in the algorithm outout although. Colleting \(M\) samples of a signal, which CoRoT did for us, is equivalent to cutting a finite piece of the real signal, i.e., we take a window between \(T_1\) and \(T_2\) of a defined signal of \({-\infty}\) to \({+\infty}\). There is a possibility that this window perfectly represents the oscillations of the original signal, that is contain complete periods of all frequencies, thus we will have a FFT without distortions, but it is statistically unlikely that this will happen. In this case, we will have incomplete periods, which will introduce high frequency distorcions in the signal that, in general, will be found at the ends of the signal.
The reason that distorcions are commonly introduced at the ends of the signal is justified by the way FFT works. Their recursive calculations imply append the end of the signal with it’s beginning (from the FFT’s premise, we would have a complete peridod), and it is at this conection that unexpected frequencies are introduced. So, to reduce this effect, we will apply the method of expanding the signal edges. Given a sequence \(A\), of size \(k\), the method consists of adding \(N\) points before \(A\)’s first point, \(A[0]\) and \(N\) points after its last point, \(A[-1]\). Finnaly, we will have a sequence \(B\), of size \((k + 2N)\).
[1]:
def expand_borders(array, num):
aux_pre = np.zeros(num)
aux_pos = np.zeros(num)
i = 0
for i in range(num):
aux_pre[i] = array[0]
aux_pos[i] = array[-1]
return np.concatenate((aux_pre, array, aux_pos)).ravel()
Avoiding wraparound error
The product \([H(u)F(u)]\), provided by the basis filtering equation, is obtained by the Array product. Consider the following 4x1 vectors:
The Array product is defined by:
In mathematics, convolution is a linear operator of two functions (\(f\) and \(g\)) that measures how the shape of one is modified by the other. The convolution of two functions involves rotating by 180° one about its origin and displacemening it past the other. In discrete functions, it is given by:
where \(f\) and \(g\) are arrays of size \(n\).
The 1-D Convolution theorem, in the case of discrete convolution, tells us that:
where \(\mathcal{F}\) is the Fourier transform and \(k\) is a constant.
Therefore, the Fourier transform of the convolution of two functions in spatial domain corresponds to the product of the transformations of these functions in the frequency domain. Writing in order to facilitate the relationship between ther Convolution theorem and the basic filtering equation, we have:
The double arrow means that the expression on the right side is obtained by taking the Fourier transform of the expression on the left, as well as, the expression on the left was obtained by taking the inverse Fourier transform of the expression on the right.
Assim, se utilizar a DTF e o Teorema da convolução para obter o mesmo resultado de filtragem, devemos tomar conta a periodicidade implícita na expressão para a DFT. Nesse caso, estaríamos aplicando a convolução à duas funções teóricamente períodicas, o que implicaria no resultado da convolução ser também periódico que, por meio de aplicações práticas nota-se ser completamente errôneo. A proximidade dos períodos de ambas funções é tamanha que eles interferem um no outro, provocando o erro de wraparound (efeito de borda).
Já estamos minimizando esse erro realizando a expansão das extremidades do sinal, contudo essa prática pode não ser suficiente para lidar com o efeito de borda. Para tornar nosso algoritmo ainda mais robusto quanto ao erro de wraparound e otimizar o cálculo da FFT, aplicaremos o Zero Padding
Esse procedimento consiste em adicionar zeros no final do sinal do dominío espacial, antes de ser calculada sua transformada de Fourier, visando aumentar o seu tamanho. Esse procedimento é utilizado para remover a periodicidade implícita do sinal, o que acarreta em resultados inesperados no resultado da filtragem. Tal propriedade matemática está deduzida na obra Fast Fourier Transform and Its Applications by E. Brigham (1988).
Outra grande utilidade do Zero Padding é em otimizar o algoritmo da FFT. Como, seguindo algumas regras, podemos modificar o tamanho no nosso sinal original, é uma boa prática transformar o seu tamanho em uma potência de 2. Para prosseguirmos com esse procedimento, temos que seguir algumas regras: FFT Zero Padding.
[2]:
def padding(array):
return np.append(array, np.zeros(len(array)))
Centering the transform
With Zero Padding applied, the next step is to try to simplift the definition of the filter transfer function, \(H(u)\). Based on symmetry properties of funcions, it is noted that the specification of \(H(u)\) is quite simplified when we are working with symmetric functions by their center, which requires that their Fourier transform \(F(u)\) also be centralized.
The centralization is demostraded by the propert of translation and rotation of the Fourier transform:
That is, multiply \(f(x)\) by the exponential shown shifts the origin of DFT to \(u_0\). Taking \(u_0 = M/2\), the exponential becomes \(e^{j\pi x} = (-1)^{x}\), stricly because \(x\) is an integer. So,
Therefore, the multiplication of the input signal \(f(x)\) by \((-1)^{x}\) implies the displacement of the data so that \(F(0)\) it in the center of the interval \([0, M-1]\), as we wanted.
This procedure must be done before computing the Fourier transform of \(f(x)\) and it is useful to help in visualizing the filtering process and in optimizing the generation of the filter array, so it not mandatory to implement it.
[3]:
def multiplying_by_minus_one_to_index(array):
i = 0
new_array = np.ones(len(array))
for i in range(len(array)):
new_array[i] = array[i] * ( (-1)**(i) )
return new_array
Taking the Fourier transform
[4]:
def fourier_transform(array):
fft = np.fft.fft(array)
return fft
Filter Transfer Function
Each filter has their own transfer function, …
Inverse Fourier transform
[5]:
def inverse_fourier_transform(array):
ifft = np.fft.ifft(array)
return ifft
Removing Zero-Padding
[6]:
def no_padding(array):
return array[:int(len(array)/2)]
Removing borders expanded
[7]:
def remove_expanded_borders(array, param):
aux = np.delete(array, np.s_[:param])
removed = np.delete(aux, np.s_[-param:])
return removed