Filtering Filter out high frequency noise The goal here is to define a signal, compute its fourier series and then either filter the series (to eliminate high frequency noise) or compress it (to send it more efficiently). An inititial signal is defined as the function f. Other sample signals appear at the end. First define the signal. > > f:=x->exp(-x^2/10)*(sin(2*x)+2*cos(4*x)+0.4*sin(x)*sin(50*x)); f := x -> 2 exp(- 1/10 x ) (sin(2 x) + 2 cos(4 x) + .4 sin(x) sin(50 x)) Bascially, the signal f involves a smoothly varying term (given first) together with some high frequency noise (the term involving sin(50x)). Let's see how Fourier series filters out the high frequency noise. First , a plot. > plot(f,-Pi..Pi,numpoints=200,thickness=0); #Plot the signal > (1/(2*Pi))*Int(f(x),x=-Pi..Pi); a0:=evalf(%); #the coefficient of a_0 Pi / | 1/2 | | / -Pi 2 exp(- 1/10 x ) (sin(2 x) + 2 cos(4 x) + .4 sin(x) sin(50 x)) dx/Pi a0 := -.009434900925 > > a:=n->evalf( (1/Pi)*Int(f(x)*cos(n*x),x=-Pi..Pi) ); # the value of a_n Pi / | | f(x) cos(n x) dx | / -Pi a := n -> evalf(----------------------) Pi > a(1); .02293963232 > b:=n->evalf( (1/Pi)*Int(f(x)*sin(n*x),x=-Pi..Pi) ); #the value of b_n Pi / | | f(x) sin(n x) dx | / -Pi b := n -> evalf(----------------------) Pi A good exercise here is to guess which frequencies lead to the dominant Fourier coefficients. Then, the guess can be checked against the computation of a_n and b_n given in the following code. > #Evaluate some sample a(n), b(n) > for n from 10 to 12 do > n; a(n); b(n); > od; > 10 -.004924890271 -.0006500494149 11 .003717132221 .0004809832411 12 -.002918355597 -.0003662307315 > Now set up the a partial Fourier series that filters out the high frequency coefficients. > n:='n'; temp:=a0; #set up the fourier partial sum n := n temp := -.009434900925 > for n from 1 to 5 do > temp:=temp+a(n)*cos(n*x)+b(n)*sin(n*x); > od; > temp := -.009434900925 + .02293963232 cos(x) + .1307208634 sin(x) temp := -.009434900925 + .02293963232 cos(x) + .1307208634 sin(x) - .04292048328 cos(2 x) + .7540210070 sin(2 x) temp := -.009434900925 + .02293963232 cos(x) + .1307208634 sin(x) - .04292048328 cos(2 x) + .7540210070 sin(2 x) + .2814246872 cos(3 x) + .1361792985 sin(3 x) temp := -.009434900925 + .02293963232 cos(x) + .1307208634 sin(x) - .04292048328 cos(2 x) + .7540210070 sin(2 x) + .2814246872 cos(3 x) + .1361792985 sin(3 x) + 1.496270756 cos(4 x) - .01729604802 sin(4 x) temp := -.009434900925 + .02293963232 cos(x) + .1307208634 sin(x) - .04292048328 cos(2 x) + .7540210070 sin(2 x) + .2814246872 cos(3 x) + .1361792985 sin(3 x) + 1.496270756 cos(4 x) - .01729604802 sin(4 x) + .2802146281 cos(5 x) + .006936761750 sin(5 x) > n:='n'; n := n > S:=temp; #the partial fourier series of through 5th order frequency S := -.009434900925 + .02293963232 cos(x) + .1307208634 sin(x) - .04292048328 cos(2 x) + .7540210070 sin(2 x) + .2814246872 cos(3 x) + .1361792985 sin(3 x) + 1.496270756 cos(4 x) - .01729604802 sin(4 x) + .2802146281 cos(5 x) + .006936761750 sin(5 x) Now plot the partial Fourier series (stored in the variable S) and compare it with the plot of the original function. > plot(S,x=-Pi..Pi,numpoints=200,thickness=0); #plot of the 5th order fourier sum > plot(f(x),x=-Pi..Pi,numpoints=200,thickness=0); #plot of the original signal > plot({f(x),S},x=-Pi..Pi,numpoints=200,thickness=0); #plot of both together Data Compression Here the goal is to compute the Fourier series of a given signal as done above. Then we eliminate all the small Fourier coefficients (set them equal to zero) to obtain a compressed signal. By changing the tolerance, we can change the degree to which the signal has been compressed. Then we can plot the compressed signal and compare it to the original signal. > > First define the signal h > h:=x->exp(-x^2/10)*(sin(2*x)+2*cos(4*x)+0.4*sin(x)*sin(10*x)); h := x -> 2 exp(- 1/10 x ) (sin(2 x) + 2 cos(4 x) + .4 sin(x) sin(10 x)) > plot(h,-Pi..Pi); #plot the signal > (1/(2*Pi))*Int(h(x),x=-Pi..Pi); a0:=evalf(%); #the coefficient of a_0 Pi / | 1/2 | | / -Pi 2 exp(- 1/10 x ) (sin(2 x) + 2 cos(4 x) + .4 sin(x) sin(10 x)) dx/Pi a0 := -.009374277785 > a:=n->evalf( (1/Pi)*Int(h(x)*cos(n*x),x=-Pi..Pi) ); # the value of a_n Pi / | | h(x) cos(n x) dx | / -Pi a := n -> evalf(----------------------) Pi > b:=n->evalf( (1/Pi)*Int(h(x)*sin(n*x),x=-Pi..Pi) ); #the value of b_n Pi / | | h(x) sin(n x) dx | / -Pi b := n -> evalf(----------------------) Pi > n:='n'; #free up the value of n and then store the Fourier coefficients n := n > for n from 1 to 25 do > aa(n):=a(n); bb(n):=b(n); > od: > > Now we ignore all the Fourier coefficients that are small (smaller than some tolerance that is specified in the variable tol) by setting them equal to zero. We store the new coefficients in the arrays aa(n) and bb(n) > n:='n'; tol:=.05; # free up n and the value of tolerance set to .05 n := n tol := .05 > for n from 1 to 25 do > if abs(aa(n)) < tol then aa(n):=0 fi; > if abs(bb(n)) < tol then bb(n):=0 fi; > od: > #Print the Fourier Coefficients > for n from 1 to 25 do > n; aa(n); bb(n); > od; 1 0 .1307208634 2 0 .7540210070 3 .2812153477 .1361792985 4 1.496587431 0 5 .2796707116 0 6 0 0 7 0 0 8 0 0 9 .1606226668 0 10 0 0 11 -.1500244133 0 12 0 0 13 0 0 14 0 0 15 0 0 16 0 0 17 0 0 18 0 0 19 0 0 20 0 0 21 0 0 22 0 0 23 0 0 24 0 0 25 0 0 > > > n:='n'; # free up the value of n n := n > hcond:=a0+sum(aa(n)*cos(n*x)+bb(n)*sin(n*x),n=1..25); # the condensed Fourier series hcond := -.009374277785 + .7540210070 sin(2 x) + 1.496587431 cos(4 x) + .1307208634 sin(x) + .1606226668 cos(9 x) + .2812153477 cos(3 x) - .1500244133 cos(11 x) + .1361792985 sin(3 x) + .2796707116 cos(5 x) > plot({h(x),hcond},x=-Pi..Pi, numpoints=200); #Plot the condensed fourier series Fast Fourier Transform The next two sections also filter and compress a given signal, but this time the Fast Fourier Transform is used to compute the Fourier coefficients. Filter out high frequency noise First define the signal. > f:=x->exp(-x^2/10)*(sin(2*x)+2*cos(4*x)+0.4*sin(x)*sin(50*x)); f := x -> 2 exp(- 1/10 x ) (sin(2 x) + 2 cos(4 x) + .4 sin(x) sin(50 x)) > plot(f,0..2*Pi, thickness=0); > readlib(FFT); #load FFT package proc(m, x, y) ... end > with(linalg): #load linear algebra package Warning, new definition for norm Warning, new definition for trace A value of N must be chosen for the FFT, which corresponds to a discretization of the signal into 2^N data points. The step size is also chosen (to be the interval length divided by 2^N). > N:=8; step:=2*Pi/2^N; N := 8 step := 1/128 Pi Now discretize the function into 2^N data points and assign to an array called func. The array ifunc is also defined which contains the imaginary parts of the input signal. Of course, most signals are real valued and in these cases, the values of ifunc are all set to zero. > func:=array([seq(evalf(f((j*step))),j=1..2^N)]): > ifunc:=array([seq(0,j=1..2^N)]): #imaginary values=0 > FFT(N,func,ifunc): #compute FFT with 2^N data points The FFT computes the Discrete Fourier Transform and stores the real and imaginary parts of these coefficients in the arrays func and ifunc (the original values are lost). > n:='n'; n := n > for n from 250 to 256 do #Now examine some of them > n; func[n]; ifunc[n]; > od; 250 .0985468415 17.51892531 251 .0098290607 25.56351769 252 10.87658358 55.41761976 253 105.2793012 -5.428184162 254 -16.99080599 -40.56820728 255 5.939679009 41.96121501 256 31.65406937 -1.856771874 The high frequency coefficients are the ones centered about n=2^N/2 (i.e. the center of the range of coefficients). To filter out high frequency, set all Fourier coefficients above say, n=7 and below 2^N - 7 equal to zero. Of course, the value 7 is arbitrary and can be changed depending on what is considered high frequency. > for n from 7 to 249 do > func[n]:=0; ifunc[n]:=0; > od: > Now compute the inverse FFT using the package iFFT. The real and imaginary parts of the signal values will be computed and stored in the arrays func and ifunc (the Fourier coefficients will now be lost). > iFFT(N,func,ifunc): > To plot this filtered signal we arrange a sequence of [x,y] values with x representing the x-coordinate of the graph (so x=j*step) and y represents the y - coordinate of the graph (so y=func[j]). > filtergraph:=seq([evalf(j*step),func[j]],j=1..2^N): > plot([filtergraph],style=line, thickness=0); Now compare with the plot of the original signal > plot(f,0..2*Pi,numpoints=500,thickness=0); > n:='n'; n := n > Data Compression Now we want to compress a given signal by computing its FFT and setting all FFT coefficients that are small equal to zero. Then we'll recompute the signal using the inverse FFT and compare with the original. > > N; step; #check N and step size 8 1/128 Pi First define the signal f > f:=x->exp(-x^2/10)*(sin(2*x)+2*cos(4*x)+0.4*sin(x)*sin(10*x)); f := x -> 2 exp(- 1/10 x ) (sin(2 x) + 2 cos(4 x) + .4 sin(x) sin(10 x)) > plot(f,0..2*Pi, thickness=0); > step; #check the step size 1/128 Pi As before, discretize the function using step size interval/2^N > func:=array([seq(evalf(f(j*step)),j=1..2^N)]): > ifunc:=array([seq(0,j=1..2^N)]): #Imaginary parts =0 > readlib(FFT); #load FFT package proc(m, x, y) ... end > FFT(N,func, ifunc): #compute FFT > > > j:='j'; j := j The tolerance variable (tol) is assigned a tolerance and a variable count is initialized to zero. The variable count will ultimately contain the number of Fourier coefficients which are set to zero. > tol:=4; count:=0; tol := 4 count := 0 Now set all Fourier coefficients which are less than tol=2 equal to zero > for j from 1 to 2^N do > if abs(func[j]) < tol then func[j]:=0; count:=count+1; fi; > if abs(ifunc[j])< tol then ifunc[j]:=0; count:=count+1; fi; > od: > > j:='j'; count; j := j 459 Note that the number of nonzero coefficients (both real and imaginary) which are NON-zero is 2(2^N)-count. Examine some Fourier coefficient values > for j from 1 to 10 do > j; func[j]; ifunc[j]; > od; 1 20.22875183 0 2 31.65241412 0 3 5.935172353 -41.88972558 4 -17.00098762 40.68825381 5 105.2586326 5.618055447 6 10.83585862 -55.11387875 7 0 -25.04167212 8 0 -16.44758304 9 0 -9.737370813 10 11.98235898 -10.89739009 > > iFFT(N,func,ifunc); #Compute inverse FFT 256 > > Set up the graphs for the function computed by condensing the FFT values and for the original. > condgraph:=seq([evalf(j*step),func[j]],j=1..2^N): > > plot([condgraph],style=line, thickness=0); > plot(f,0..2*Pi,style=line, thickness=0); > > > > > noise:=x->2*x*exp(-(x^2)); p:=x->-52*x^4+100*x^3-49*x^2+2; 2 noise := x -> 2 x exp(-x ) 4 3 2 p := x -> -52 x + 100 x - 49 x + 2 > f:=x->p(x)+noise((x-1/3)*100)+noise((x-2/3)*200); f := x -> p(x) + noise(100 x - 100/3) + noise(200 x - 400/3) > with(plots): > p2:=plot(f,0..1,numpoints=500, thickness=0): > display({p2}, thickness=0); > > Here is another example of a signal that can be used for filtering and data compression. It consists of a polynomial p defined over the inteval 0 to 1 with some noise added at x=1/3 and x=2/3. Since the interval is now [0,1] care must be taken to change the step size and plot intervals if these functions are used in the above. > noise:=x->2*x*exp(-(x^2)); p:=x->-52*x^4+100*x^3-49*x^2+2; > f:=x->p(x)+noise((x-1/3)*100)+noise((x-2/3)*200); > > > p:=x->a*x^4+b*x^3+c*x^2+2; 4 3 2 p := x -> a x + b x + c x + 2 > eq1:=p(1/2)=-1; eq1 := 1/16 a + 1/8 b + 1/4 c + 2 = -1 > eq2:=D(p)(1/2)=0; eq2 := 1/2 a + 3/4 b + c = 0 > eq3:=p(1)=1; eq3 := a + b + c + 2 = 1 > solve({eq1,eq2,eq3},{a,b,c}); {c = -49, b = 100, a = -52} > assign(%); > plot(p,0..1); > p(x); 4 3 2 -52 x + 100 x - 49 x + 2 > plotsetup(x11); >