# Filtering using fourier series

Suppose I have a measured signal $$M$$ that has frequency components from 0 to 50 Hz. I plot the specturm of this signal using FFT and I observe its frequency content (power vs frequency). Then, I decompose signal $$M$$ to its fourier series components and weights, I call this series $$M_{hat}$$. I compare the spectrum of $$M_{hat}$$ and $$M$$ and they look identical.

Now, I decide to extract the the frequency content between 0 to 10 Hz from $$M_{hat}$$. I consider the sins and cosines from 0 to 10Hz linearly combined with their weights. I call this signal $$M_{hat10}$$. When I plot the spectrum of $$M_{hat10}$$ and $$M_{hat}$$ on top of each other, I don't see them exactly matching in the range 0 to 10Hz. However, $$M_{hat10}$$ is indeed cut after 10Hz.

I don't understand why don't $$M_{hat}$$ and $$M_{hat10}$$ match in the range 0 to 10Hz.

As per OPs real signal shared in txt file, it looks like extraction was done only one 1 side corresponding to FFT $$M_{hat}$$. For a real signal, we need to do extraction from other side also. So effectively we need to combine weights of double the number of components. The spectrum of the new signal $$M_{hat10}$$ will match that of $$M_{hat}$$ in the frequency range 0 to 10Hz as shown below

clc
clear all
close all

A=table2array(T);
n=A(:,1);  %time
x=A(:,2);  %signal

M=fft(x);
L=length(M);
Z=round(10/50*L); % computing FFT points corresponding to 10Hz
%Extracting only Z points from either side of FFT and setting rest to zero
Mhat10 = [M(1:Z); zeros(L-2*Z,1); M(L-Z+1:L)];

xhat10 = ifft(Mhat10,L);
Mhat10_fft=fft(xhat10,L);

figure()
plot(abs(M))
hold on
plot(abs(Mhat10_fft))


• My original signal (M) length is 3000 and my approximated signal "Mhat" length is also 3000. I have decomposed Mhat recursively giving a matrix of sins and cosines with different frequencies each with length 3000 --> I got the weights. When I wanted to "cut" from Mhat to get Mhat10, I multiplied the corresponding sines and cosines with the weights I got. My sines and cosines have length 3000. – HaneenSu Apr 12 '20 at 6:58
• And how many Fourier series components do you compute? – jithin Apr 12 '20 at 7:01
• I computed 2500 components for the sins and another 2500 components for the cosins – HaneenSu Apr 12 '20 at 7:03
• I tried by making signal length and sin/cosine length different and they are not matching. But when they are same it is matching for frequencies 0-10Hz. Have you tried that? – jithin Apr 12 '20 at 7:22
• I implemented your code exactly on my signal but it doesn't seem to work. My data is from this link:motusbioengineering.com/testdata/chfig8.txt . The first column is time and the second column is the signal. – HaneenSu Apr 12 '20 at 7:51

This will match exactly if the time domain signal is properly recreated from the Inverse DFT as expressed in the sine and cosine form. In this case the DFT and inverse DFT are transform pairs so can be expressed in either domain with no loss of information.

Here is the general form for (painfully) computing the inverse DFT from sines and cosines which helps show the mathematical utility of always using the complex form of frequency when dealing with complex sinals ($$e^{j\omega t}$$ instead of $$\cos(\omega t)+j\sin(\omega t)$$). (And ultimately using the FFT algorithm and not the direct DFT equation is much simpler but this is insightful).

Given the Inverse DFT as:

$$x_n = \sum_{k=0}^{N-1}X_ke^{j2\pi nk/N}$$

When N is odd this is equally of the form:

$$x_n = X_0 + \sum_{k=1}^{(N-1)/2}\big(X_ke^{j2\pi nk/N}+ X_{[N-k]}e^{-j2\pi nk/N}\big)$$

If $$x_n$$ was real then $$X_{[N-k]} = X^*_k$$ leading to further simplifications, but in general for complex $$x_n$$ the above is given as below using $$\omega_n = 2\pi n/N$$ and the relationship $$A e^{j\theta} = A\cos(\theta)+jA\sin(\theta)$$:

$$x_n = X_0+\sum_{k=1}^{(N-1)/2}(X_k+X_{[N-k]})cos(k \omega_n)+j\sum_{k=1}^{(N-1)/2}(X_k-X_{[N-k]})sin(k\omega_n)$$

When N is even the inverse DFT is of the form:

$$x_n = X_0+\sum_{k=1}^{N/2-1}\big(X_ke^{j2\pi nk/N}+ X_{[N-k]}e^{-j2\pi nk/N}\big)+ X_{[N/2]}e^{j\pi n}$$

Which similarly becomes:

$$x_n = X_0 +(-1)^nX_{[N/2]} \ldots \\+ \sum_{k=1}^{N/2-1}(X_k+X_{[N-k]})cos(k \omega_n)+j\sum_{k=1}^{N/2-1}(X_k-X_{[N-k]})sin(k\omega_n)$$

In the OP's case where only the first $$m$$ frequencies of the DFT are used (corrresponding to $$X_0$$ when $$m=0$$, $$X_1$$ and $$X_{[N-1]}$$ when $$m=1$$ etc... and notably $$X_{[N/2]}=0$$, then regardless of $$N$$ is odd or even the result is:

$$x_n = X_0 + \sum_{k=1}^m(X_k+X_{[N-k]})cos(k \omega_n)+j\sum_{k=1}^m(X_k-X_{[N-k]})sin(k\omega_n)$$

I also want to add this important additional note that this is actually a very poor way to implement a filter (it is called the frequency sampling approach, where the coefficients of the implemented FIR filter would be the IFFT of this frequency domain mask). It will indeed provide the results exactly as given in frequency at each bin center (in this case pass each frequency at the bin centers below 10 Hz with magnitude 1, and each frequency above 10 Hz with magnitude 0), but for all frequencies in between the bin centers both in the passband and stopband it will have a lot more deviation from other optimized filter design approaches (such as least squares and equi-ripple filter design).

This response at this post by @hotpaw2 explains this further in more detail: Why is it a bad idea to filter by zeroing out FFT bins?.

This is because once you consider only the linear combinations of sins and cosines in the 0 to 10Hz, you have essentially taken away from the energy and representation of the original time domain signal in the 0 to 50Hz. In essence this new signal formed from the sins and cosines from 0 to 10Hz essentially represents a different time domain signal itslef. Hence the fourier contents will not match. The convergence between $$M_{hat10}$$ and $$M_{hat}$$ in 0 to 10Hz will only be in the mean square sense.