DSP Project 2 Report
FFT Related Studies BY Hani Saleh, Adel Hussein and Bassam Mohd
1
Table of Contents Table of Contents ............................................................................................................... 2 Table of Figures ................................................................................................................. 4 List of Tables ...................................................................................................................... 4 1
Abstract....................................................................................................................... 5
2
Introduction................................................................................................................ 5
3 FFT Leakage, Uncertainty Principle, Filtering Using FFT & Arithmetic Fourier transform ............................................................................................................................ 6 3.1
Technical background.................................................................................................6
3.1.1 The Discrete Fourier Transform ............................................................................................. 6 3.1.2 The Discrete Fourier Transform Spectral Leakage................................................................. 6 3.1.3 Using Weighted Windows to Minimize Leakage ................................................................... 9 The term window has a specific meaning when applied to samples used to create a Fourier spectrum. Consider the signal in the first figure. In order to analyze the signal for harmonics we have to sample it for some fixed period (T)................................................................................................................... 9 3.1.3.1 Hanning Window Example.............................................................................................. 10 3.1.4 Time-Frequency Uncertainty Principle................................................................................. 13 3.1.5 Arithmetic Fourier Transform .............................................................................................. 13 3.1.6 Filtering using Fast Fourier Transform................................................................................. 16
3.2
Design detailed info ................................................................................................16
3.2.1 3.2.2
3.3
Arithmetic Fourier Transform .............................................................................................. 16 Filtering using Fast Fourier Transform................................................................................. 17
Computer simulation..............................................................................................17
3.3.1 3.3.2 3.3.3 3.3.4
FFT Leakage......................................................................................................................... 17 FFT filtering ......................................................................................................................... 21 AFT....................................................................................................................................... 24 Uncertainty ........................................................................................................................... 26
3.4
References...............................................................................................................27
3.5
Hani Saleh Matlab Code ........................................................................................28
3.5.1 3.5.2 3.5.3 3.5.4 3.5.5 3.5.6 3.5.7 3.5.8 3.5.9 3.5.10
4
Arithmetic Fourier Transform .............................................................................................. 28 Test the AFT......................................................................................................................... 29 Mobius Function................................................................................................................... 31 Test the Mobius function ...................................................................................................... 32 Is Square Free function......................................................................................................... 34 Test the FFT Leakage ........................................................................................................... 35 The FFT filter kernel function .............................................................................................. 39 Test the FFT filtering............................................................................................................ 41 Test Uncertainty ................................................................................................................... 49 RADIX-2 FFT.................................................................................................................. 51
Radix4 FFT (Adel Hussein Module)....................................................................... 54 4.1
Introduction:...........................................................................................................54
4.2
Background: ...........................................................................................................55
4.3
Design detailed information: ..................................................................................63
4.4
Computer simulation..............................................................................................65
2
5
4.5
Performance Evaluation.........................................................................................70
4.6
Conclusion ..............................................................................................................71
4.7
References...............................................................................................................72
4.8
Matlab Code ...........................................................................................................73
Fast Fourier Transform: Matrix Form (Bassam Jamil Module) .......................... 79 5.1
FFT Matrix-Form Theory......................................................................................79
5.2
Computation of FFT Matrix-Form........................................................................81
5.2.1
A’s and D’s Matrices............................................................................................................ 83
5.3
MatLab program Flow...........................................................................................84
5.4
Results .....................................................................................................................87
5.4.1 5.4.2 5.4.3 5.4.4
5.5 5.5.1 5.5.2 5.5.3 5.5.4
4-point Example from Professor Sos Notes.......................................................................... 87 8-point Example from Professor Sos Notes.......................................................................... 87 Compare Matrix Form FFT with Matlab Functions ............................................................. 88 Number of operations ........................................................................................................... 91
Matlab Code ...........................................................................................................91 Bmohd_fft_matrix_top.m ..................................................................................................... 92 bmohd_fft_matrix_form ....................................................................................................... 94 bmohd_A_kn.m .................................................................................................................... 96 bmohd_calc_operations ........................................................................................................ 96
3
Table of Figures Figure 1: Period extension in DFT...................................................................................... 7 Figure 2: Sinusoid signal spectrum (no Leakage)............................................................... 8 Figure 3: Complex signal plot (No windowing) ................................................................. 9 Figure 4: Complex signal plot (Samplind window)............................................................ 9 Figure 5: Hanning window example................................................................................. 11 Figure 6: Time window shapes ......................................................................................... 12 Figure 7: AFT Flowchart .................................................................................................. 15 Figure 8: Blackman window............................................................................................. 18 Figure 9: Hamming window ............................................................................................. 18 Figure 10: Hanning window ............................................................................................. 19 Figure 11: Blackman window........................................................................................... 20 Figure 12: All windows combined on the same figure ..................................................... 21 Figure 13: FFT filtering, low pass .................................................................................... 22 Figure 14: FFT filtering, High pass .................................................................................. 22 Figure 15: FFT filtering, Band Pass.................................................................................. 23 Figure 16: FFT filtering, Band stop .................................................................................. 23 Figure 17: FFT filtering, all filters shown together .......................................................... 24 Figure 18: AFT vs FFT spectrums.................................................................................... 25 Figure 19: AFT reconstructed Time signal....................................................................... 26 Figure 20: Uncertainty ...................................................................................................... 27 Figure 21 Layers of Matrix Computation ......................................................................... 82 Figure 22 Permutation Matrix (N=8)................................................................................ 83 Figure 23 A's and D's Matrices ......................................................................................... 84 Figure 24 Program Flow ................................................................................................... 87 Figure 25 FFT: Matrix Form vs. Matlab........................................................................... 91 List of Tables Table 1, Window functions comparision .......................................................................... 12 Table 2, Error measurements for the FFT filters .............................................................. 24 Table 3, Error measurements for the AFT algorithm ....................................................... 26
4
1 Abstract The main objective of this project is to implement 3 fast algorithms for the DFT: Radix4, Matrix form and the Arithmetic Fourier Transform (AFT). Also we studied the effect of FFT leakage, Time - Frequency Uncertainty relationship and Low pass, High Pass, Band Pass and Band stop filtering using the FFT.
2 Introduction Chapters 3, 4 & 5 list the results of the work we have done on this project. In the following chapters each individual module will be presented. • Ch3: Hani Saleh module o Arithmetic Fourier Transform o FFT Leakage o Filtering Using FFT o Time-Frequency Uncertainty relationship • Ch4: Adel Hussein module: o FFT Radix4 Implementation • Ch5: Bassam Mohd module: o FFT implementation in Matrix form
5
3 FFT Leakage, Uncertainty Principle, Filtering Using FFT & Arithmetic Fourier transform (Hani Saleh Module) These modules have been developed by Hani Saleh. They include: o DFT Leakage o Uncertainty Principle o Filtering Using FFT o Arithmetic Fourier Transform 3.1 Technical background The following technical background is based on the references attached and a numerous amount of data collected form the World Wide Web. 3.1.1
The Discrete Fourier Transform
The Fourier transform is a time domain to frequency domain mathematical transform. For any given continuous signal the Fourier transform provides a decomposition of the signal into the amplitudes and frequencies of a set of sine waves which would reproduce the original signal when summed. This property makes it easy to identify the frequencies at which most of the signal strength is being transmitted, making the Fourier transform an ideal tool for spectrum analysis. To perform the transform by machine the signal has to sampled and the discrete Fourier transform (DFT) used. For an N sample transform the DFT is defined by the formula:
Where x(n) is the sample at time index of n and X(k) is a vector of N values at frequency index k corresponding to the magnitude of the sine waves at resulting from the decomposition of the time indexed signal. 3.1.2
The Discrete Fourier Transform Spectral Leakage
When you use the DFT/FFT to find the frequency content of a signal, it is inherently assumed that the data that you have is a single period of a periodically repeating waveform. This is shown in the following figure. The first period shown is the one sampled. The waveform corresponding to this period is then repeated in time to produce the periodic waveform.
6
Figure 1: Period extension in DFT [Ref-7] As seen in the previous figure, because of the assumption of periodicity of the waveform, discontinuities between successive periods will occur. This happens when you sample a non integer number of cycles such as 7.5. These artificial discontinuities turn up as very high frequencies in the spectrum of the signal, frequencies that were not present in the original signal. These frequencies could be much higher than the Nyquist frequency, and will be aliased somewhere between 0 and fs/2. The spectrum you get by using an FFT therefore will not be the actual spectrum of the original signal, but will be a smeared version. It appears as if the energy at one frequency has leaked out into all the other frequencies. This phenomenon is known as spectral leakage.
The left hand side of the following figure shows a sine wave sampled. The timedomain waveform has an integer number of cycles, so the assumption of periodicity does not create any discontinuities.
7
Figure 2: Sinusoid signal spectrum (no Leakage)
In the right hand side in the figure above, you see the spectral representation when you sample a non integer number of cycles of the time waveform. The periodic extension of this signal creates a discontinuity similar to periodic extension figure. Notice how the energy is now spread over a wide range of frequencies, so the relative height difference between the FFT peak amplitude and the neighboring bins is reduced. This smearing of the energy is spectral leakage. Leakage is a direct result of using a limited set of sampled data. If we could collect and process an infinite amount of samples, the frequency increment (xaxis) becomes very small and all possible frequencies would be valid harmonics. The FFT algorithm does not add anymore leakage than other algorithms but it makes correcting the problem more difficult by limiting our choices for N to powers of 2. There are two approaches to improving the accuracy of the spectrum when leakage is a problem. 1. Change the sample set to line up with the harmonic content. 2. Modify the sampled values with a weighted window.
8
3.1.3
Using Weighted Windows to Minimize Leakage
The term window has a specific meaning when applied to samples used to create a Fourier spectrum. Consider the signal in the first figure. In order to analyze the signal for harmonics we have to sample it for some fixed period (T).
Figure 3: Complex signal plot (No windowing) [Ref-8] Mathematically, this amount to multiplying one set of values by 1.0 and all other values of the waveform by zero. In effect, we create a set of sampled data by blocking out all other parts of the waveform with a Rectangular Window with an amplitude 1.0. The result is shown in the next figure.
Figure 4: Complex signal plot (Samplind window) [Ref-8]
9
Most computer algorithms, including the FFT, process this set of samples as if it were a continuous loop. If N=256 (samples 0 to 255), then it assumes that the value after 255 is the same as the value at sample 0. This assumption is true when the window (T) includes an integer multiple of a repeating cycle. When it is not true, as in this example, the discontinuity between the value of the first and last samples looks like a step change to the FFT and causes leakage. Using a non-rectangular window can reduce leakage. In practice, this is done after the sample set is collected by scaling each sample in the set. The scaled values are than used as the input to the FFT. Many windows have been developed and each has its own set of costs and benefits. 3.1.3.1 Hanning Window Example Figure 4 shows the effect of applying a Hanning window to a pure sine tone. The left-top plot shows a sine tone that is not periodic in the time window without the windowing function resulting in leakage in the FFT (left-bottom). When a Hanning window is applied (top-right), then the leakage is reduced in the FFT (bottom-right). The resulting spectrum is a sharp narrow peak with amplitude of 1.0. Notice that it does not have exactly the same shape as the FFT of the original periodic sine wave in Figure 3, but the amplitude and frequency errors resulting from leakage are corrected. A Windowing function minimizes the effect of leakage to better represent the frequency spectrum of the data.
10
Figure 5: Hanning window example [Ref-9]
Notes on Selecting a Window 1. If the input signal is periodic and a multiple of the sampling window, then use a Rectangular window. 2. If the input is a short pulse or burst that starts and ends at the same amplitude, then use a Rectangular window (as long as the sampling window includes the entire transient). 3. If the input is a sections of a continuous waveform that is not periodic then use the Hamming or Triangular windows.
The most common windows and their features are given below. This table can be used to choose the best windowing function for each application.
11
Table 1, Window functions comparison [Ref-9]
Figure 6: Time window shapes [Ref-9]
12
3.1.4 Time-Frequency Uncertainty Principle There is an inverse relationship between the dispersion of a function and the range of the frequencies which are present in its transform. Thus the shorter is the duration of a transient signal, the wider is the spread of the frequencies in its transform. In mathematical physics, an analogous relationship between the spatial dispersion of a wave train and its frequency dispersion is the basis of the uncertainty principle of Heisenberg. Recall that Fourier basis elements b (t) = eit exhibit poor time localization abilities - a consequence of the fact that b (t) is evenly spread over all t 2 (−1,1). By time localization we mean the ability to clearly identify signal events which manifest during a short time interval, such as the "glitch" described in an earlier example1. At the opposite extreme, a basis composed of shifted Dirac deltas b(t) = delta(t −T) would have excellent time localization but terrible "frequency localization," since every Dirac basis element is evenly spread over all Fourier frequencies . By frequency localization we mean the ability to clearly identify signal components which are concentrated at particular Fourier frequencies, such as sinusoids. These observations motivate the question: Is there a basis that provides both excellent frequency localization and excellent time localization? The answer is "not really": there is a fundamental tradeoff between the time localization and frequency localization of any basis element. The desire to have highly localized basis functions to handle transients, with localized Fourier transforms to obtain good coding gain, are contradictory requirements due to the Heisenberg uncertainty relation between a function and its Fourier transform. The selection of the basis functions must be a compromise between these conflicting requirements. 3.1.5 Arithmetic Fourier Transform [Ref-1] Assume that we want to represent the real valued periodic signal A(t) by Fourier series:
13
Ak & θk are the amplitude and the phase of the Kth harmonic respectively. Let us consider a bank of N delay-line (or tranversal) filters each with input a(t). The output of the nth filter is denoted by S(n,t) and defined by the formula:
Note that S(1,t) = A(t). By replacing the variable t by (t+1/n) and using periodicity of A(t), it can be seen that S(n,t) is periodic with with period 1/n dor n=1,2,…., N. Further, by calculating the frequency response of the linear filtering operation of the above formula we can write the Fourier series of S(n,t) as:
Because of our assumption that all of the harmonics above N are zero, then only the first [N/n] terms of the above formula can be non zero. [N/n] denotes the largest integer that is less than or equal to N/n. It can be proven that [Ref-1] the Fourier coefficients can be given by:
Where µ(m) is the Mobius (1904) function defined as: µ1(1) = 1, µ1(n) = (-1)s if n = (p1) (p2) (p3)..(pS), where pi are distinct primes, If n is a square-free positive integer with an even or odd number of distinct prime factors. A number is said to be squarefree (or sometimes quadratfree; Shanks 1993) if its prime decomposition contains no repeated factors. 14
µ1(n) = 0
if p2|n for any prime p (n is not square-free)
The following figure shows a flowchart for computing the AFT:
Figure 7: AFT Flowchart [Ref-1]
15
3.1.6 Filtering using Fast Fourier Transform The DFT generate the spectrum of the input signal; the frequency content of the input signal will be present in the signal. Digital filtering can be easily done with the FFT spectrum, by: 1) Calculate the FFT of the input signal. 2) Chang the FFT results values which correspond to the frequencies you wish to change 3) Perform IFFT to get the filtered time domain signal.
3.2
Design detailed info
3.2.1 Arithmetic Fourier Transform The following algorithm was coded in Matlab code to calculate the AFT using the formulas given in the above theory section: 1) Calculate the moving average values and stored in an array bn_0(1…N/2), N2 = N/2: for n= 1:N2 for m = 1:(1*n -0) indx1 = N -(N*(m)/(1*n) ) + 1; indx2 = (N/(2*n)); indx1_norm = floor(indx1); interp = indx1 - indx1_norm; if (indx1_norm > N) indx1_norm = N; end if (indx2_norm > N) indx2_norm = N; end if((interp ==0)|| (indx1_norm == N)) bn_0(n) = bn_0(n) + (1)^m*din(indx1_norm); else if(interp > 0) bn_0(n) = bn_0(n) + (1)^m*(din(indx1_norm) + din(indx1_norm +1))/2; else bn_0(n) = bn_0(n) + (1)^m*(din(indx1_norm) + din(indx1_norm -1))/2; end end end
16
bn_0(n) = bn_0(n) / (1 * n); end 2) Calculate the Fourier coefficients data; a(1…n): for n = 1:N2 for k = 1:(1*floor(N2/n)) a(n) = a(n) + mobius(k)*bn_0(1*k*n); end y(n+1) = N2*a(n); y(N-n+1) = N2*a(n); end 3.2.2 Filtering using Fast Fourier Transform Since the FFT generate the spectrum of the input signal; it was used to perform digital filtering by getting ride off the unneeded frequency components. Digital filtering using FFT was performed as follows: • Calculate the FFT of the input signal. • Generate a frequency response in the frequency domain that corresponds to the Filter requirements • Multiply the filter vector by the FFT result of the input signal • Perform IFFT to get the filtered time domain signal. 3.3
Computer simulation
3.3.1 FFT Leakage To study the leakage effect a spectrum of a sinusoid has been generated without applying any windows and then the spectrum of the same signal has ben generated using Blackman, Hamming and Hanning windows. Results are shown in the following figures:
17
Figure 8: Blackman window
Figure 9: Hamming window
18
Figure 10: Hanning window
The following figure shows the plot of time domain and the spectrum for the Blackman, Hamming and the Hanning windows.
19
Figure 11: Blackman window
The following figure shows all of the windows reduction in leakage side by side:
20
Figure 12: All windows combined on the same figure Conclusion: Applying windows makes the input signal symmetric and results with much less leakage. The Hamming & Hanning windows reduction in leakage was better than the Blakman. 3.3.2 FFT filtering Matlab code has been developed to: 1) Create the signal s1= cos(2*pi*25*t1) + cos(2*pi*200*t1) +cos(2*pi*400*t1); 2) Performed a. Low pass filtering with cutoff frequency = 100 Hz b. High pass filtering with cutoff frequency = 100 Hz c. Band pass filtering with pass band frequency = 100 - 300 Hz d. Band stop filtering with cutoff frequency = 100 - 300 Hz The annotated plots of the input signals and the filtered output are shown in the following figures:
21
Figure 13: FFT filtering, low pass
Figure 14: FFT filtering, High pass
22
Figure 15: FFT filtering, Band Pass
Figure 16: FFT filtering, Band stop
23
Figure 17: FFT filtering, all filters shown together The following table summarizes the error measurements for the above FFT filters:
Error-Max Error-RMS SNR
Error Measurements for the Above FFT filters Low High BP 0.835961 0.835601 0.789634 0.042446 0.0422272 0.0446035 24.4513 27.4956 24.01
BS 0.789634 0.0446035 27.0252
Table 2, Error measurements for the FFT filters Conclusion: Filtering using FFT was very effective in achiving the desired output. Due to the need of calculating the FFT of the input signal and then the need to calculate the IFFT of the filtered spectrum such an approach might not be suitable for high speed real tim signal processing applications.. 3.3.3 AFT The FFT algorithm has been implemented using Matlab. The spectrum for the signal: S1 = cos(2*pi*5*t1)+cos(2*pi*10*t1)+cos(2*pi*15*t1);
24
Has been generated using the AFT algorithm and the Matlab FFT algorithm. The following figure shows the time domain signal and the spectrum of both the AFT algorithm and the FFT algorithm.
Figure 18: AFT vs FFT spectrums The following figure shows the original time domain signal and the AFT reconstructed time signal.
25
Figure 19: AFT reconstructed Time signal
Error-Max Error-RMS SNR
FFT 4.44089e-016 2.34776e-016 314.391
AFT , 0.971114 0.287952 12.6182
Table 3, Error measurements for the AFT algorithm Conclusion: The AFT algorithm used showed some promising results with fairly high error margin. The main reason for the error is the use of zero interpolation to calculate the values of the missing samples due to the non uniform sampling requirement of the AFT algorithm. 3.3.4 Uncertainty To study the uncertainty a Gaussian function has been generated using Matlab. The width of the time domain has been varied and the spectrum has been studied. The following figure shows the time domain signal and the resultant spectrum.
26
Figure 20: Uncertainty Conclusion: There is an inverse relationship between the dispersion of a function in the time domain and the range of the frequencies which are present in its Fourier transform. The simulation results confirmed that the shorter is the duration of a signal in time domain, the wider is the spread of the frequencies in its Fourier transform. 3.4 1.
2. 3. 4. 5. 6. 7.
References “The arithmetic Fourier transform”, Sadasiv, G.; ASSP Magazine, IEEE [see also IEEE Signal Processing Magazine] , Volume: 5 , Issue: 1 , Jan. 1988 Pages:13 - 17H. Nyquist, "Certain topics in telegraph transmission theory,"
Trans. AIEE, vol. 47, pp. 617-644, Apr. 1928. Professor SOS, DSP Class notes, University of Texas at San Antonio, Fall 2004. Oppenheim, Schafer, Buck, Discrete-Time Signal Processing, Second Edition, Prentice Hall 1999. “Inside The Fft Black Box- Serial And Parallel Fast Fourier Transform Algorithms”, Crc Press - Chu And George Et Al Digital Signal Processing, Schaum’s Outlines Understanding Digital Signal Processing, 1st edition, Richard Lynos http://zone.ni.com/devzone/conceptd.nsf/webmain/4D2F2B3C170CC02E862568 C1006857EB
27
8. 9.
3.5 3.5.1
http://www.astro-med.com/knowledge/fourier.html http://www.ldsgroup.com/docs/site_documents/AN014_Understanding_FFT_2.pdf
Hani Saleh Matlab Code Arithmetic Fourier Transform
function d=aft(din) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % aft_5.m % Radix-2 FFT, Decimation in Frequency % % Custom Arithmetic Fourier transform.' % aft_5(din) is the discrete Fourier transform (DFT) of vector data_in. % %This function has been developed by Hani Saleh % For DSP project 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% no_add = 0; no_mult = 0; N = length(din); N2 = floor(N/2); M = log2(N); % remove the dc component from the input signal a0 = mean (din); din = din - a0; %Initialize the data a(1:N)=0; bn_0(1:N) = 0; y(1:N) = 0; % Calculate the moving average values for n= 1:N2 for m = 1:(1*n -0) indx1 = N -(N*(m)/(1*n) ) + 1; indx2 = (N/(2*n)); indx1_norm = floor(indx1); interp = indx1 - indx1_norm; indx2_norm = floor(indx1+indx2+0.5); if (indx1_norm > N) indx1_norm = N; end if (indx2_norm > N) indx2_norm = N; end if((interp ==0)|| (indx1_norm == N)) bn_0(n) = bn_0(n) + (1)^m*din(indx1_norm); else
28
if(interp > 0) bn_0(n) = bn_0(n) + (1)^m*(din(indx1_norm) + din(indx1_norm +1))/2; else bn_0(n) = bn_0(n) + (1)^m*(din(indx1_norm) + din(indx1_norm -1))/2; end end end bn_0(n) = bn_0(n) / (1 * n); end % Create the Fourier output real vector for n = 1:N2 for k = 1:(1*floor(N2/n)) a(n) = a(n) + mobius(k)*bn_0(1*k*n); end y(n+1) = N2*a(n); y(N-n+1) = N2*a(n); end %Return the results d= y; % %Print adds 7 Mult no report % cmplx=N*log(N); % no_add; % no_mult; % fprintf('-----------------\n'); % fprintf('N^2, NLog(N) , No of Adds , No of Mult\n'); % fprintf('%g , %g , %g , %g\n',N^2,cmplx,no_add,no_mult); % fprintf('-----------------\n');
3.5.2
Test the AFT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % This m-file test the Arithmetic Fourier Transfrom % % % Written by Hani Saleh for DSP project 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Creat 4 time axisbases, at four sampling rates. fs1 = 100; % Define the number of samples num_samples = 100; %Create the time sampling vectors t1 = 0:(1/fs1):num_samples*(1/fs1)-(1/fs1); % Four samples of a signal that is composed of two cosines functions
29
freq1 = 5; freq2 = 10; freq3 = 15; s1 = cos(2*pi*freq1*t1)+cos(2*pi*freq2*t1)+cos(2*pi*freq3*t1); f_axis1 = (-1/2)*fs1:(fs1/num_samples):(1/2)*fs1-(fs1/num_samples); s1_bar = s1 - mean(s1); tmr1 = cputime; a1 = aft(s1); tmr2 = cputime - tmr1 tmr1 = cputime; a2 = fft(s1); tmr3 = cputime - tmr1 S1 = fftshift(abs(a1)); S2 = fftshift(abs(a2)); aft_rec = ifft(a1); fft_rec = ifft(a2); diff = s1 - aft_rec; diff1 = s1 - fft_rec; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % %%%%%%%%%%%%%%%%% s1 error measurment s1_error = s1 - aft_rec; s1_error_max = max (abs(s1_error)); s1_error_rms = std (s1_error); s1_snr = 10*log10 ( (std(s1))^2/ (std(s1_error))^2); %%%%%%%%%%%%%%%%% s2 error measurment s2_error = s1 - fft_rec; s2_error_max = max (abs(s2_error)); s2_error_rms = std (s2_error); s2_snr = 10*log10 ( (std(s1))^2/ (std(s2_error))^2); out_s1 = sprintf(' ERROR Measurements, Error = Original-Signal IFFT-of-AFT/FFT results \n'); out_s1a = sprintf(' AFT \t\t FFT \n'); out_s2 = sprintf('Error-Max = %g \t %g \t \n', s1_error_max,s2_error_max); out_s3 = sprintf('Error-RMS = %g \t %g \t \n', s1_error_rms,s2_error_rms); out_s4 = sprintf('SNR = %g \t %g \t \n', s1_snr,s2_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s1a, out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % plot the time domain signals plot(aft_rec,'g'); hold on; plot(s1,'r');
30
title('Input Signal and the AFT Reconstructed signal'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Plots of the four samples figure; subplot(2,2,1); plot(t1,s1); %axis([0 0.1 -3 3]); grid; title('Input Signal'); subplot(2,2,2); plot(t1,s1_bar); %axis([0 0.1 -3 3]); grid; title('Input Signal - mean value'); subplot(2,2,3); plot(f_axis1,S1); %axis([0 0.1 -3 3]); grid; title('AFT Generated Spectrum'); subplot(2,2,4); plot(f_axis1,S2); %axis([0 0.1 -3 3]); grid; title('Matlab FFT Generated Spectrum'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3.5.3 Mobius Function function dout=mobius(din) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % mobius function % mobius(n)= u(n) is defined for all positive natural numbers n and has its values in {-1, 0, 1} % depending on the factorization of n into prime factors. It is defined as follows % % u(n) = 1 if n is a square-free positive integer with % an even number of distinct prime factors. % u(n) = -1 if n is a square-free positive integer with % an odd number of distinct prime factors. % u(n) = 0 if n is not square-free. % This is taken to imply that u(1) = 1. The value of u(0) is generally left undefined, % but the Maple computer algebra system for example returns ?1 for this % values so this function returns -1 for 0 as well. %%%%%%%%%%%%%%%%%% % This function returns the mobius function result for % an input integer %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31
% Test if the input is an integer >=0 is_int = (abs(round(din)) - abs(din)); if ( is_int ~= 0) fprintf('-----------------\n'); fprintf('Error, the input number is not an integer %g \n',din); fprintf('-----------------\n'); dout= NaN; return elseif (din < 0) fprintf('-----------------\n'); fprintf('Error, the input number is a negative integer %g \n',din); fprintf('-----------------\n'); dout= NaN; return end %Check if it is square free number or not fact = factor(din); leng = length(fact); sq_free = is_sq_free(din); %Return the results if(din == 1) dout= 1; elseif (din == 0) dout = -1; elseif (sq_free == 1) dout = (-1)^leng; else dout= 0; end 3.5.4 Test the Mobius function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%% % This script tests the mobius function % developed by Hani saleh for the DSP % project 2 %%%%%%%%%%%%% % the following 3 sequences were taken from % http://www.maa.org/editorial/mathgames/mathgames_11_03_03.html % Which gives the results of the mobius function for certain integers % % This program compares the results returned by the function to the % expected resulted and flags an error if it does not match
32
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%mobius = 0 not_sq_free= [4,8,9,12,16,18,20,24,25,27,28,32,36,40,44,45,48,49,50,52,54,56,60,63,64,68,72,75,76,8 0,81,84,88,90,92,96,98,99,100,104,108,112,116,117,120,121,124,125,126,128,132,135,1 36,140,144,147,148,150,152,153,156,160]; %mobius = -1 sq_free_odd=[2,3,5,7,11,13,17,19,23,29,30,31,37,41,42,43,47,53,59,61,66,67,70,71,73,7 8,79,83,89,97,101,102,103,105,107,109,110,113,114,127,130,131,137,138,139,149,151,1 54,157,163,165,167,170,173,174,179,181,182,186,190,191,193]; %mobius = 1 sq_free_even=[1,6,10,14,15,21,22,26,33,34,35,38,39,46,51,55,57,58,62,65,69,74,77,82,8 5,86,87,91,93,94,95,106,111,115,118,119,122,123,129,133,134,141,142,143,145,146,155 ,158,159,161,166,177,178,183,185,187,194,201,202,203,205,206,209,210,213,214]; no_errors = 0; for i = 1:length(not_sq_free) if (mobius(not_sq_free(i)) ~= 0) fprintf('-----------------\n'); fprintf('Error number - %g - is not square free and result is 0\n',not_sq_free(i)); fprintf('-----------------\n'); no_errors = no_errors + 1; end end for i = 1:length(sq_free_odd) if (mobius(sq_free_odd(i)) ~= -1) fprintf('-----------------\n'); fprintf('Error number - %g - is odd square free and result is not -1\n',sq_free_odd(i)); fprintf('-----------------\n'); no_errors = no_errors + 1; end end for i = 1:length(sq_free_even) if (mobius(sq_free_even(i)) ~= 1) fprintf('-----------------\n'); fprintf('Error number - %g - is odd square free and result is not +1\n',sq_free_even(i)); fprintf('-----------------\n'); no_errors = no_errors + 1;
33
end end if (no_errors ~= 0) fprintf('-----------------\n'); fprintf('The total number of errors = %g \n',no_errors); fprintf('-----------------\n'); else fprintf('-----------------\n'); fprintf('Congratulations your function is working \n',no_errors); fprintf('-----------------\n'); end 3.5.5 Is Square Free function function dout=is_sq_free(din) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % square free number % A number is said to be squarefree (or sometimes quadratfrei; Shanks 1993) % if its prime decomposition contains no repeated factors. All primes are therefore % trivially squarefree. The number 1 is by convention taken to be squarefree. % The squarefree numbers are 1, 2, 3, 5, 6, 7, 10, 11, 13, 14, 15, ... (Sloane's A005117). % The squareful numbers (i.e., those that contain at least one square) are 4, 8, 9, 12, 16, 18, % 20, 24, 25, ... (Sloane's A013929). %%%%%%%%%%%%%%%%%% %This function returns one if the passed number is square free % and 0 otherwise %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %din =15; fact = factor(din); leng = length(fact); sq_free = 0; if (leng >= 2) for i = 1:leng fact_rpt = find(fact == fact(i)); no_rptns = length(fact_rpt); if (no_rptns > 1) sq_free = 0;
34
break; else sq_free = 1; end end else sq_free = 1; end %sq_free; %Return the results dout= sq_free; 3.5.6
Test the FFT Leakage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % This m-file illustrates the effects of DFT leackage % and windows enhancments to FFT results % in time and frequency domains % % % Written by Hani Saleh for DSP project 2 % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fprintf('Study the effect of the sampling frequency on the freq. domain\n'); fprintf('of the sampled signal.\n'); % Creat 4 time axisbases, at four sampling rates. fs1 = 5; % Define the number of samples num_samples = 32; %Create the time sampling vectors t1 = 0:(1/fs1):num_samples*(1/fs1)-(1/fs1); % Four samples of a signal that is composed of two cosines functions freq1 = 1; win1 = blackman(num_samples)'; win2 = hamming(num_samples)'; win3 = hann(num_samples)'; s1 = cos(2*pi*freq1*t1) ; s2 = s1 .* win1; s3 = s1 .* win2; s4 = s1 .* win2;
35
f_axis1 = (-1/2)*fs1:(fs1/num_samples):(1/2)*fs1-(fs1/num_samples); S1 S2 S3 S4
= = = =
fftshift(abs(fft_radix2_hsaleh(s1,num_samples))); fftshift(abs(fft_radix2_hsaleh(s2,num_samples))); fftshift(abs(fft_radix2_hsaleh(s3,num_samples))); fftshift(abs(fft_radix2_hsaleh(s4,num_samples)));
WIN1 = fftshift(abs(fft_radix2_hsaleh(win1,num_samples))); WIN2 = fftshift(abs(fft_radix2_hsaleh(win2,num_samples))); WIN3 = fftshift(abs(fft_radix2_hsaleh(win3,num_samples))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% PLOTTING Blackman %%%%%%%%%%%%%%%%%%%%%%%% % Plots of the four samples figure; subplot(2,2,1); plot(t1,s1); %axis([0 0.1 -3 3]); grid; title('Input signal time domain; Freq=1,FS=5'); %%%%%%% subplot(2,2,2); plot(t1,s2); %axis([0 0.1 -3 3]); grid; title('Input * blackman window'); %%%%%%%%%%%%%% % Plots of the four windows subplot(2,2,3); plot(f_axis1,S1,'-x'); hold on; stem(f_axis1,S1,'-x'); %axis([0 0.1 -3 3]); grid; title('Input spectrum no window'); %%%%% subplot(2,2,4); plot(f_axis1,S2,'-x'); hold on stem(f_axis1,S2,'-x'); %axis([0 0.1 -3 3]); grid; title('Input spectrum blackman window'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% PLOTTING Hamming %%%%%%%%%%%%%%%%%%%%%%%% % Plots of the four samples figure; subplot(2,2,1); plot(t1,s1); %axis([0 0.1 -3 3]); grid; title('Input signal time domain; Freq=1,FS=5'); %%%%%%% subplot(2,2,2);
36
plot(t1,s3); %axis([0 0.1 -3 3]); grid; title('Input * Hamming window'); %%%%%%%%%%%%%% % Plots of the four windows subplot(2,2,3); plot(f_axis1,S1,'-x'); hold on; stem(f_axis1,S1,'-x'); %axis([0 0.1 -3 3]); grid; title('Input spectrum no window'); %%%%% subplot(2,2,4); plot(f_axis1,S3,'-x'); hold on stem(f_axis1,S3,'-x'); %axis([0 0.1 -3 3]); grid; title('Input spectrum Hamming window'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% PLOTTING Hanning %%%%%%%%%%%%%%%%%%%%%%%% % Plots of the four samples figure; subplot(2,2,1); plot(t1,s1); %axis([0 0.1 -3 3]); grid; title('Input signal time domain; Freq=1,FS=5'); %%%%%%% subplot(2,2,2); plot(t1,s3); %axis([0 0.1 -3 3]); grid; title('Input * Hanning window'); %%%%%%%%%%%%%% % Plots of the four windows subplot(2,2,3); plot(f_axis1,S1,'-x'); hold on; stem(f_axis1,S1,'-x'); %axis([0 0.1 -3 3]); grid; title('Input spectrum no window'); %%%%% subplot(2,2,4); plot(f_axis1,S4,'-x'); hold on stem(f_axis1,S4,'-x'); %axis([0 0.1 -3 3]); grid; title('Input spectrum Hanning window'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% PLOTTING all together %%%%%%%%%%%%%%%%%%%%%%%% % Plots of the four samples figure; subplot(2,4,1); plot(t1,s1); %axis([0 0.1 -3 3]); grid; title('Input signal time domain; Freq=1,FS=5'); %%%%%%% subplot(2,4,2); plot(t1,s2); %axis([0 0.1 -3 3]); grid; title('Input * blackman window'); %%%%%%%% subplot(2,4,3); plot(t1,s3); %axis([0 0.1 -3 3]); grid; title('Input * hamming window'); %%%%%%%% subplot(2,4,4); plot(t1,s4); %axis([0 0.1 -3 3]); grid; title('Input * hanning window'); %%%%%%%%%%%%%% % Plots of the four windows subplot(2,4,5); stem(f_axis1,S1,'-x'); %axis([0 0.1 -3 3]); grid; title('Input spectrum no window'); %%%%% subplot(2,4,6); stem(f_axis1,S2,'-x'); %axis([0 0.1 -3 3]); grid; title('Input spectrum blackman window'); %%%%%% subplot(2,4,7); stem(f_axis1,S3,'-x'); %axis([0 0.1 -3 3]); grid; title('Input spectrum hamming window'); %%%%%%% subplot(2,4,8); stem(f_axis1,S4,'-x'); %axis([0 0.1 -3 3]); grid; title('Input spectrum hanning window'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
38
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Plots of windows %%%%%%%%%%%%%%%%%%% figure; subplot(2,3,1); plot(t1,win1); %axis([0 0.1 -3 3]); grid; title('Blackman window Time domain'); %%%% subplot(2,3,2); plot(t1,win2); %axis([0 0.1 -3 3]); grid; title('Hamming window Time domain'); %%%%%% subplot(2,3,3); plot(t1,win3); %axis([0 0.1 -3 3]); grid; title('Hanning window Time domain'); % Plots of the four windows %%%%%% subplot(2,3,4); plot(f_axis1,WIN1); axis([-2 2 0 15]); grid; title('Blackman window spectrum'); %%%%%% subplot(2,3,5); plot(f_axis1,WIN2); axis([-2 2 0 20]); grid; title('Hamming window spectrum'); %%%%%%% subplot(2,3,6); plot(f_axis1,WIN3); axis([-2 2 0 20]); grid; title('Hanning window spectrum');
3.5.7 The FFT filter kernel function function d=filter_kern_fft(data,f_sample, f_cutoff1,f_cutoff2,type) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% % filt_kern_fft.m % Filtering using FFT % Written by Hani Saleh for the DSP project 2 %----------------------------------------------------------------% Input: % data = data to be filtered % fsample = sampling frequency % f_cutoff1 = cutoff frequency for low pass and high pass and it is % the first edge of the band pass and band reject 39
% f_cutoff2 = The second edge of the band pass and band reject, % set to 0 for low pass and high pass % % type = type of filter % lo = low pass % ho = high pass % bp = bandpass % br = band reject %---------------------------------------------------------------%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Written by Hani Saleh % For the DSP class project 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % define variables N = length(data); fs = f_sample; fc1 = f_cutoff1; fc2 = f_cutoff2; N2 = int16(N/2); fc1_num = int16(N * (fc1/fs)); fc2_num = int16(N * (fc2/fs));
% low pass if (type == 'lo') data(fc1_num:N2)=0; data(N2+1:N-fc1_num)=0; end % high pass if (type == 'hi') data(1:fc1_num)=0; data(N-fc1_num:N)=0; end % bandpass if (type == 'bp') data(1:fc1_num-1)=0; data(fc2_num+1:N2-1)=0; data(N2:N-fc2_num-1)=0; data(N-fc1_num+1:N)=0; end
40
% bandpass if (type == 'bs') data(fc1_num:fc2_num)=0; data(N-fc2_num:N-fc1_num)=0; end
%Return the results d= data; 3.5.8 Test the FFT filtering %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % This m-file illustrates the filtering using FFT % % % Written by Hani Saleh for DSP project 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Creat 4 time axisbases, at four sampling rates. fs1 = 1000; % Define the number of samples num_samples = 1024; %Create the time sampling vectors t1 = 0:(1/fs1):num_samples*(1/fs1)-(1/fs1); %High sampling rate time vector for the time domain plot t_plot = 0:0.001:0.1; % Four samples of a signal that is composed of two cosines functions freq1 = 25; freq2 = 200; freq3= 400; signal1 = cos(2*pi*freq1*t1); signal2 = cos(2*pi*freq2*t1); signal3 = cos(2*pi*freq3*t1); s1 = signal1 + signal2 + signal3; s_plot = cos(2*pi*freq1*t_plot) + 2*cos(2*pi*freq2*t_plot);
S1 = fft_radix2_hsaleh(s1,num_samples);
41
f_axis1 = (-1/2)*fs1:(fs1/num_samples):(1/2)*fs1-(fs1/num_samples); %filter the data S1_filt = filter_kern_fft(S1,fs1,100,0,'lo'); S2_filt = filter_kern_fft(S1,fs1,100,0,'hi'); S3_filt = filter_kern_fft(S1,fs1,100,300,'bp'); S4_filt = filter_kern_fft(S1,fs1,100,300,'bs'); % reconstruct the time domain data s1_reconst = ifft_radix2_hsaleh(S1_filt,num_samples); s2_reconst = ifft_radix2_hsaleh(S2_filt,num_samples); s3_reconst = ifft_radix2_hsaleh(S3_filt,num_samples); s4_reconst = ifft_radix2_hsaleh(S4_filt,num_samples); s1_orig = signal1; s2_orig = signal2+signal3; s3_orig = signal2; s4_orig = signal1+signal3; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%% s1 error measurment s1_error = s1_orig - s1_reconst; s1_error_max = max (abs(s1_error)); s1_error_rms = std (s1_error); s1_snr = 10*log10 ( (std(s1_orig))^2/ (std(s1_error))^2); %%%%%%%%%%%%%%%%% s2 error measurment s2_error = s2_orig - s2_reconst; s2_error_max = max (abs(s2_error)); s2_error_rms = std (s2_error); s2_snr = 10*log10 ( (std(s2_orig))^2/ (std(s2_error))^2); %%%%%%%%%%%%%%%%% s3 error measurment s3_error = s3_orig - s3_reconst; s3_error_max = max (abs(s3_error)); s3_error_rms = std (s3_error); s3_snr = 10*log10 ( (std(s3_orig))^2/ (std(s3_error))^2); %%%%%%%%%%%%%%%%% s4 error measurment s4_error = s4_orig - s4_reconst; s4_error_max = max (abs(s4_error)); s4_error_rms = std (s4_error); s4_snr = 10*log10 ( (std(s4_orig))^2/ (std(s4_error))^2); out_s1 = sprintf(' ERROR Measurements, Error = Original-Signal - IFFT-of-Radix2Form-FFT \n'); out_s1a = sprintf(' Low \t\t High \t\t BP \t\t BS\n');
42
out_s2 = sprintf('Error-Max = %g \t %g \t %g \t %g \t \n', s1_error_max,s2_error_max,s3_error_max,s4_error_max); out_s3 = sprintf('Error-RMS = %g \t %g \t %g \t %g \t \n', s1_error_rms,s2_error_rms,s3_error_rms,s4_error_rms); out_s4 = sprintf('SNR = %g \t %g \t %g \t %g \t \n', s1_snr,s2_snr,s3_snr,s4_snr); Errror_measures = sprintf('%s %s %s %s',out_s1,out_s1a, out_s2, out_s3, out_s4); display (Errror_measures); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots for the low pass figure; %%%%%%%%%%%%%% subplot(2,2,1); % stem(t1,s1,'-*'); % hold on; plot(t1,s1,'-g'); % hold off; axis([0 0.1 -3 3]); grid; title('Input signal'); %%%%%%%%%%%%%% subplot(2,2,3); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S1)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input Signal Spectrum FFT'); %%%%%%%%%%%%%% subplot(2,2,2); % stem(t4,s4); % hold on; plot(t1,abs(s1_reconst),'-r'); hold on; plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('Low pass Reconstructed signal (green- Original, Red - Reconst'); %%%%%%%%%%%%%% subplot(2,2,4); %stem(t3,s3);
43
% hold on; plot(f_axis1,fftshift(abs(S1_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Low pass Filtered signal fft');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots for the high pass figure; %%%%%%%%%%%%%% subplot(2,2,1); % stem(t1,s1,'-*'); % hold on; plot(t1,s1,'-g'); % hold off; axis([0 0.1 -3 3]); grid; title('Input signal'); %%%%%%%%%%%%%% subplot(2,2,3); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S1)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% subplot(2,2,2); % stem(t4,s4); % hold on; plot(t1,abs(s2_reconst),'-r'); hold on; plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('High pass Reconstructed signal'); %%%%%%%%%%%%%% subplot(2,2,4);
44
%stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S2_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('High pass Filtered signal fft'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots for the BandPass pass figure; %%%%%%%%%%%%%% subplot(2,2,1); % stem(t1,s1,'-*'); % hold on; plot(t1,s1,'-g'); % hold off; axis([0 0.1 -3 3]); grid; title('Input signal'); %%%%%%%%%%%%%% subplot(2,2,3); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S1)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% subplot(2,2,2); % stem(t4,s4); % hold on; plot(t1,abs(s3_reconst),'-r'); hold on; plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('Band pass Reconstructed signal'); %%%%%%%%%%%%%% subplot(2,2,4);
45
%stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S3_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Band pass Filtered signal fft'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots for the Band Stop pass figure; %%%%%%%%%%%%%% subplot(2,2,1); % stem(t1,s1,'-*'); % hold on; plot(t1,s1,'-g'); % hold off; axis([0 0.1 -3 3]); grid; title('Input signal'); %%%%%%%%%%%%%% subplot(2,2,3); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S1)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% subplot(2,2,2); % stem(t4,s4); % hold on; plot(t1,abs(s4_reconst),'-r'); hold on; plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('Band Stop Reconstructed signal'); %%%%%%%%%%%%%% subplot(2,2,4);
46
%stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S4_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Band Stop Filtered signal fft'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % % Plots of all the samples combined figure; %%%%%%%%%%%%%% subplot(2,4,1); % stem(t1,s1,'-*'); % hold on; plot(t1,signal1,'-g'); % hold off; axis([0 0.1 -1.5 1.5]); grid; title('Signal1, Freq = 25 Hz'); %%%%%%%%%%%%%% subplot(2,4,4); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S1)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% subplot(2,4,2); % stem(t4,s4); % hold on; plot(t1,signal2,'-r'); % hold on; % plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('Signal2, Freq = 200 Hz'); %%%%%%%%%%%%%% subplot(2,4,3);
47
% stem(t4,s4); % hold on; plot(t1,signal3,'-r'); % hold on; % plot(t1,abs(s1),'-g'); axis([0 0.1 -1.5 1.5]); hold off; % axis([0 0.1 -3 3]); grid; title('Signal3, Freq = 400 Hz'); %%%%%%%%%%%%%% subplot(2,4,5); %stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S1_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Low pass Filtered signal fft'); %%%%%%%%%%%%%% subplot(2,4,6); %stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S2_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('High pass Filtered signal fft'); %%%%%%%%%%%%%% subplot(2,4,7); %stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S3_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Band pass Filtered signal fft'); %%%%%%%%%%%%%% subplot(2,4,8); %stem(t3,s3); % hold on; plot(f_axis1,fftshift(abs(S4_filt)),'-g'); % hold off; % axis([0 0.1 -3 3]);
48
grid; title('Band reject Filtered signal fft'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 3.5.9
Test Uncertainty
%sampling freq fs1 = 100; % Define the number of samples num_samples = 256; %Create the time sampling vectors t1 = 0:(1/fs1):num_samples*(1/fs1)-(1/fs1); %%%%%%%%%%%%%%%%%%%%%%%%% % y = gaussmf(x,[sig c]) %%%%%%%%%%%%%%%%%%%%%%%%% c= 0.5; sigma1 = 1; sigma2 = 0.1; sigma3 = 0.01; sigma4 = 0.001; z1 = (-.5*(t1 - c).^2)/(sigma1^2); z2 = (-.5*(t1 - c).^2)/(sigma2^2); z3 = (-.5*(t1 - c).^2)/(sigma3^2); z4 = (-.5*(t1 - c).^2)/(sigma4^2); y1=(1/(sqrt(2*pi)))*exp(z1); y2=(1/(sqrt(2*pi)))*exp(z2); y3=(1/(sqrt(2*pi)))*exp(z3); y4=(1/(sqrt(2*pi)))*exp(z4); f_axis1 = (-1/2)*fs1:(fs1/num_samples):(1/2)*fs1-(fs1/num_samples); S1=fft(y1); S2=fft(y2); S3=fft(y3); S4=fft(y4); % % Plots of the four samples figure; %%%%%%%%%%%%%% subplot(2,4,1); % stem(t1,s1,'-*');
49
% hold on; plot(t1,y1,'-g'); % hold off; % axis([-0.5 1.5 -1 2]); grid; title('Input signal'); %%%%%%%%%%%%%% subplot(2,4,5); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S1)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% %%%%%%%%%%%%%% subplot(2,4,2); % stem(t1,s1,'-*'); % hold on; %axis([-0.5 1.5 -1 2]); plot(t1,y2,'-g'); % hold off; %axis([0 0.1 -3 3]); grid; title('Input signal'); %%%%%%%%%%%%%% subplot(2,4,6); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S2)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% %%%%%%%%%%%%%% subplot(2,4,3); % stem(t1,s1,'-*'); % hold on; %axis([-0.5 1.5 -1 2]); plot(t1,y3,'-g'); % hold off; %axis([0 0.1 -3 3]); grid; title('Input signal');
50
%%%%%%%%%%%%%% subplot(2,4,7); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S3)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% %%%%%%%%%%%%%% subplot(2,4,4); % stem(t1,s1,'-*'); % hold on; %axis([-0.5 1.5 -1 2]); plot(t1,y4,'-g'); % hold off; %axis([0 0.1 -3 3]); grid; title('Input signal'); %%%%%%%%%%%%%% subplot(2,4,8); % stem(t2,s2); % hold on; plot(f_axis1,fftshift(abs(S4)),'-g'); % hold off; % axis([0 0.1 -3 3]); grid; title('Input FFT'); %%%%%%%%%%%%%% 3.5.10 RADIX-2 FFT function d=fft_radix2_hsaleh(data_io,number_samples) % fft_radix2.m % Radix-2 FFT, Decimation in Frequency % % Custom Discrete Fourier transform.' % FFT(data_io) is the discrete Fourier transform (DFT) of vector data_io. % % For length N input vector data_io, the DFT is a length N vector data_io,' % with elements' % N' % data_io(k) = sum data_io(n)*edata_iop(-j*2*pi*(k-1)*(n1)/N), 1 <= k <= N.' % n=1' % The inverse DFT (computed by IFFT) is given by' % N'
51
% data_io(n) = (1/N) sum 1)/N), 1 <= n <= N.' % k=1'
data_io(k)*edata_iop( j*2*pi*(k-1)*(n-
%Detrmine if your input is real data vector if (mean(imag(data_io))== 0) real_data = 1; else real_data = 0; end N = number_samples; M = log2(number_samples); % Compute the FFT using decimation in Freq for L = 1: M L; LE = 2^(M+1-L); LE2 = LE/2; % Wiggle factors for 1st level U = 1.0; %Compute the Wigle factors S = cos(double(pi/LE2)) - j* sin(double(pi/LE2)); for J = 1:LE2 J; for I = J:LE:N I; IP = I + LE2; T = data_io(I) + data_io(IP); data_io(IP) = (data_io(I) - data_io(IP)) * U; data_io(I) = T; end U = U * S; end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Perform the decimation in frequency ND2 = N/2; NM1 = N -1; %start with an index of 1 J = 1; for I=1:NM1 if I < J I; % Swap data_io(J) with data_io(I) T = data_io(J); data_io(J) = data_io(I); data_io(I) = T; end K = ND2; while (K < J) J = J -K; K = K /2; end
52
%Increment the index J = J + K; end %For real data input make sure than data_io(N/2+m) = data_io(m) for m=2 ... N/2 % Due to the limitation of floating point accuracy there might be very % small diffrences, this part of the function remove these if (real_data==1) for INDEdata_io = 2:N/2 data_io(N+2-INDEdata_io)=data_io(INDEdata_io)'; end end %Return the results d= data_io;
53
4 Radix4 FFT (Adel Hussein Module) The objective of project 2 is to implement different Fast Fourier Transform (FFT) algorithms to solve the Discrete Fourier Transform (DFT). The over all project will have more components to it like using windowing and the leakage effect, another will implement Matrix Form of Fast Fourier Transform (FFT). This portion of the project will discuss the implementation of the Radix-2 FFT and the Radix-4 FFT. These two algorithms that will be covered in this part. Results that the algorithms are consistent with each other will be verified by plotting the magnitude and phase spectrum of the signal for each algorithm. We will run the 8 points example in the class to show the correctness of the implementation.
4.1
Introduction:
In real-life applications the analog signal needs to be processed properly to take advantage of modern days technology. We have to develop techniques to process these signals using today’s available technology. In 1807 Fourier discovered that a wide class of signals can be made by summing scaled sine and cosine functions. This provided good tool for analyzing both periodic and stationary signals. The direct implementation of the Fourier transform is relatively slow. The Fast Fourier Transform was and is being used in areas such as pattern recognition, interpolation, image reconstruction, linear filters, and others. The FFT calculation involves complex computations and time consuming. The engineering work in signal/image processing field is dependant on developing algorithms that can calculate transforms and convolution operations in today’s computers. For an aperiodic signal with finite energy : ∞
∫ x a ( t )e
Xa(F ) =
−∞ ∞
∫ X a ( F )e
xa ( t ) =
− j 2 πFt
j 2 πFt
dt
dF
- Voltage Spectrum - Inverse for recovery of signal
−∞
For a dicrete - time signal : X( ω) =
∞
∞
∑ x( n )e − jωn ⇒ X ( f ) = ∑ x( n )e − j 2πfn
n = −∞ π
x( n ) =
n = −∞
1/ 2
1 X ( ω )e jωn dω ⇒ ∫ X ( f )e j 2 πfn df 2π −∫π −1 / 2
Relationship between analog signal and discrete - time signal : X ( f ) = Fs
∞
∑ X a [( f − k )Fs ]
k = −∞
54
As mentioned above, the need for developing faster algorithms to implement the DFT. Fast Fourier Transform is an implementation of DFT. The direct implementation of FFT requires computation as follows: 1. 2N2 evaluations of trigonometric functions. 2. 4N2 real multiplications. 3. 4N(N-1) real additions. 4. Several indexing and addressing operations. 5. Total is in the order of O(N2) In short by some manipulation, this the radix-2 implementation will have calculations of the order (NlogN) computation while the radix-4 will be 60 to 80% of the time for radix2.
4.2
Background:
The Discrete Fourier Transform is a computation tool that is necessary in many digital signal processing applications. Such applications include linear filtering, correlation analysis, and spectrum analysis. The Fourier transform is a time domain to frequency domain mathematical transform. For any continuous the Fourier transform provides a decomposition of the signal into the amplitudes and frequencies of a set of sine waves which could produce the original signal when summed. This property makes it easy to identify the frequencies at which most of the signal strength is transmitted, making the Fourier transform an ideal tool for spectrum analysis. The DFT produces a unique finiteduration sequence of a discrete-time signal in the frequency domain. By having computational efficient algorithms, signals can be digitally processed faster in the frequency domain than in the time domain. Theorems and Equations •
Discrete Fourier Transform for N samples
X( k ) = =
N −1
∑ x( n )e − j 2πkn / N
k =0 N −1
∑ x( n )WNkn
k =0
, k = 0,1,2,… , N − 1
(
, where WN = e − j 2π / N
)
Where x(n) is the sample at time n and X(k) is a vector of N values at frequency index k corresponding to the magnitude of the sine waves resulting from the decimation of the
55
time index signal. The above equation suggests a calculation complexity of O(N2) is required to get the DFT. Using a fast Fourier transform will greatly reduce this complexity. The significant work produced on fast Fourier algorithms are the CooleyTukey radix-r decimation in time (DFT) and the Sandy-Tukey radix-r decimation in frequency (DIF) Fast Fourier Transform (FFT). Radix-R FFT reduce the complexity of DFT from O(N2) to O(Nlog(N)). The most popular of these are radix-2 and radix-4. Following are the explaination of the theory of the implementation of both radix-2 and raxid-4. •
Radix-2 FFT (Decimation-In-Time) Algorithm Suppose N = 2 v. M = N / 2. L = 2 Even - numbered samples : f1(n) = x( 2n) Odd - numbered samples : f 2(n) = x( 2n + 1 ) N −1
X (k ) = ∑ x(n)WNkn
, n = 0,1,…,
N −1 2
, k = 0,1,…, N − 1
n =0
X(k) can also be expressed using the odd and even parts as : N X (k ) = G1 (k ) + G2 (k ) , k = 0,1,…, − 1 2 N N X k + = G1 (k ) − G2 (k ) , k = 0,1,…, − 1 2 2 where G1 (k ) = F1 (k ) & G2 (k ) = WNk F2 (k ) and F1 (k ) =
( N / 2 ) −1
∑ f (m)W 1
km N /2
& F2 (k ) =
m =0
, k = 0,1,…,
N −1 2
( N / 2 ) −1
∑f
2
(m)W Nkm/ 2
m=0
v
For N = 2 , the radix-2 decimation-in-time will be performed v = log2N times and the computation involves: 1. (N/2)log2N complex multiplications. 2. Nlog2N complex additions.
56
Figure 1 First step in the decimation-in-time algorithm
Figure 2 Three stages in the computation of an N = 8-point DFT
57
Figure 3 Eight-point decimation-in-time FFT algorithm. (Oppenheim 9.32 p. 662)
Figure 4 butterfly computation for the decimation-in-time FFT algorithm. (Oppenheim 9.33 p. 662)
58
Figure 5 Shuffling of the data and bit reversal.
•
Radix-4 FFT (Decimation-In-Time) Algorithm Suppose N = 4 v. M = N / 4. L = 4. l , p = 0,1,2,3; m, q = 0,1,…, N / 4 − 1; n = 4m + l ; k = ( N / 4) p + q. We get 4 subsequences : x( 4n), x(4n + 1), x(4n + 2), & x(4n + 3) 3
[
]
X (k ) = X ( p, q) = ∑ W Nlq F (l , q ) W4lp
, n = 0,1,…,
N −1 4
, p = 0,1,2,3
l =0
( N / 4 ) −1
where F (l , q) =
∑ x(l, m)WNmq/ 4 m =0
, l = 0,1,2,3 , q = 0,1,2,…,
N −1 4
x(l , m) = x(4m + l ) and N X ( p, q ) = X p + q 4
For N = 2v, the radix-4 decimation-in-time will be performed v = log4N times and the computation involves:
59
1. (3N/8)log2N complex multiplications. 2. Nlog2N complex additions.
Figure 6 Basic butterfly computation in a radix-4 FFT algorithm.
Figure 7 Decimation-intime butterfly matrix. Class notes
60
Figure 8 Sixteen-point radix-4 decimation-in-time algorithm with input in normal order and output in digit reversed order.
Assumptions Made •
We assume N not to be a prime number and to be factorable by the appropriate amount (power of 2 or 4). These assumptions are not restrictive, since we can pad any sequence with zeros to ensure a factorization of the form N = LM. Since our test signal was both a power of 2 and a power of 4 (i.e., 4096 = 212 = 46), no padding was necessary for any of the calculations.
•
The processing times will be different for computers with different hardware, and so only relative performance is important. The computer which ran all of the computations has the following specs: Intel P4 with Mobile Tech 2.4 Ghz 512 MB of RAM
61
•
In order to implement the radix-2 and radix-4 algorithms, recursion will be needed for the repeated decimation-in-time.
Matlab Functions •
FFT Discrete Fourier transform. FFT(X) is the discrete Fourier transform (DFT) of vector X. For matrices, the FFT operation is applied to each column. For N-D arrays, the FFT operation operates on the first non-singleton dimension. FFT(X,N) is the N-point FFT, padded with zeros if X has less than N points and truncated if it has more. FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the dimension DIM. For length N input vector x, the DFT is a length N vector X, with elements N
X(k) = sum x(n) * exp(-j * 2 * pi * (k - 1) * (n - 1)/N), 1 <= k <= N. n =1
The inverse DFT (computed by IFFT) is given by N
x(n) = (1/N) sum X(k) * exp( j * 2 * pi * (k - 1) * (n - 1)/N), 1 <= n <= N. k =1
The relationship between the DFT and the Fourier coefficients a and b in N/2
x(n) = a0 + sum x(n) = a(k) * cos(2 * pi * k * t(n)/(N * dt)) + b(k) * sin(2 * pi * k * t(n)/(N * dt)) k =1
is a0 = X(1)/N, a(k) = 2 * real(X(k + 1))/N, b(k) = - 2 * imag(X(k + 1))/N, where x is a length N discrete signal sampled at times t with spacing dt. •
CPUTIME CPU time in seconds. CPUTIME returns the CPU time in seconds that has been used by the MATLAB process since MATLAB started. For example: t=cputime; the function; cputime-t returns the cpu time used to run the function. The return value may overflow the internal representation and wrap around.
62
4.3
Design detailed information:
In this part of the project, we will develop the radix-2 and radix-4 FFT. We will generate an input signal to test the algorithms. We will also run the 8 point FFT example given in the class notes to validate the function. Some comparison with the Matlab built in FFT functions will be performed. These functions will also be used in the main overall project. Following are the detailed steps used to implement the different algorithms of this part of the project. All were done in the Matlab environment. 1 for 0 ≤ n ≤ B−1 1. Create the periodic input signal x(n) = , where the period 0 for B ≤ n ≤ A−1 A=64 and pulse length B=A/2. Duration will be from 0 ≤ n ≤ N-1, where N=64A. 2. Take CPU time 1. Implement radix-2 decimation-in-time algorithm. Take CPU time 2. The difference in the CPU times is the total time taken for the radix-2 algorithm. Where L = M = 64. The radix-2 decimation-in-time algorithm is divided into three steps. a. Compute the N-point DFT’s by computing N/2 two-point DFT’s. ( N / 2 ) −1
∑ x (m)W
G1 (k ) =
1
km N /2
for 0 ≤ k ≤
N −1 2
for 0 ≤ k ≤
N −1 2
m =0
( N / 2 ) −1
G2 ( k ) =
∑x
2
(m)WNkm/ 2
m =0
b. Sequence continues until we get two N/2 – point DFT’s. This is expressed by H1(k) and H2(k). H 1 (k ) = G1 (k ) H 2 (k ) = G2 (k )W Nk
,0 ≤ k ≤
N −1 2
c. Compute N-point DFT’s by combining H1(k) and H2(k). ( N / 2 ) −1
X (k ) =
∑ H (k ) + H 1
2
(k )
k =0
( N / 2 ) −1 N ) = ∑ H 1 (k ) − H 2 ( k ) 2 k =0 N X = X (k ) + X (k + ) 2
X (k +
3. Plot the magnitude and phase spectra of X(ω) for -2π<ω<2π for the radix-2 algorithm.
63
4. Take CPU time 1. Implement radix-4 decimation-in-time algorithm. Take CPU time 2. The difference in the CPU times is the total time taken for the radix-4 algorithm. The radix-4 decimation-in-time algorithm is divided into three steps. a. Decimate the N-point sequence into four subsequences by selecting L = 4 and M = N/4 with the divide and conquer approach. b. The four N/4-point DFT’s are calculated. This is expressed by H1(k), H2(k), H3(k) and H4(k). c. Compute N-point DFT’s by combining Hn(k) terms according to the radix4 decimation-in-time butterfly. This yields X(4k), X(4k+1), X(4k+2), and X(4k+3). 5. Plot the magnitude and phase spectra of X(ω) for -2π<ω<2π for the radix-4 algorithm.
FLOW CHARTS FOR CREATING THE CODE Form Periodic Sequence Form x(n)
for 0≤n≤B−1 1 0 for otherB≤n≤A−1
x(n) =
A=64, B=A/2. 0 ≤ n ≤ N-1
Radix-2 Algorithm Take the radix-2 DFT. Check for base case (no recursion). Find CPU time.
Input x(n)
Decimate x(n) by 2
Compute N/2point DFTs
Plot magnitude and phase spectrum of X(w)
Combine to form X(n)
Radix-4 Algorithm Take the radix-4 DFT. Check for base case (no recursion). Find CPU time.
Input x(n)
Decimate x(n) by 4
Compute N/4point DFTs
Plot magnitude and phase spectrum of X(w)
Combine to form X(n)
Figure 9: Flow of the code
64
4.4
Computer simulation
Following are some of the output as a result of running the algorithms we built. Some will be as data and others are as plots. The code used to generate and simulate the effects on the functions we used as an input will be included in the appropriate section.
5.1 This is the output using radix-2 to run the 8 points example in the notes by professor SOS. The results match the ones in the example.
X = Columns 1 through 4 4.0000 5.0000 + 2.4142i -2.0000 + 6.0000i 5.0000 + 0.4142i Columns 5 through 8 12.0000 5.0000 - 0.4142i -2.0000 - 6.0000i 5.0000 - 2.4142i 5.2 This is the output using radix-4 to run the 8 points example in the notes by professor SOS. Since the input has only 8 elements, we have padded it with zeros at the end. The input was:
X=[4,-3,2,0,-1,-2,3,1,0,0,0,0,0,0,0,0] The results, the points in red match the ones in the example. X = Columns 1 through 4 4.0000
0.3627 + 0.0776i 5.0000 + 2.4142i 1.3286 - 3.4531i
Columns 5 through 8 -2.0000 + 6.0000i 8.0856 + 5.6179i 5.0000 + 0.4142i 6.2230 + 5.1487i Columns 9 through 12 12.0000
6.2230 - 5.1487i 5.0000 - 0.4142i 8.0856 - 5.6179i
Columns 13 through 16 -2.0000 - 6.0000i 1.3286 + 3.4531i 5.0000 - 2.4142i 0.3627 - 0.0776i
65
5.3 The output of Radix-2 of the input signal x(n)
Figure 10: Phase and Amplitude using Radix2 5.4 The output of Radix-4 of the input signal x(n)
66
Figure 11: Phase and Amplitude using Radix4 5.5 The output of filtering, windowing and leakage using Radix-4
67
68
5.6 Tabulated Error measurements
Error type
S1 - S1_ifft_Radix4
Matlab_ff t- Radix4_fft
Error_MAX
6.89412e-015
5.36504e-013
Error_RMS
3.59504e-015
5.35799e-014
SNR
299.979
300.548
5.7 The Radix-4 compared to Matlab FFT functiom
69
Figure 12: Using Radix4 on input
4.5
Performance Evaluation
Each method of the Fourier Transform was checked for correct magnitude and phase spectra, and the resulting plots are given. More importantly, the results may be summarized by various methods to show time improvements due to the various FFT algorithms Several key phenomena are apparent in the tabulated and graphed results. • Each Fourier Transform method is progressively faster, with the original mathematical definition being far slower than any improved algorithm. • The fastest method is Matlab’s FFT algorithm, also having the smallest time to completion of any of the methods. • Each method reduces the number of complex multiplications and/or additions over the previous method. 1.
The decimation-in-time Radix-2 FFT algorithm (Radix-2). •
The magnitude and phase spectra are as expected.
•
An improvement over the Direct and Divide and Concor. The recursive nature of this method is the key to the large performance gain.
•
The Radix-2 FFT algorithm should be implemented with a calculation-in-place system, improving the speed further. In our case, we followed the theory closely in order to achieve the most straightforward computation.
2.
The decimation-in-time Radix-4 FFT algorithm (Radix-4). •
The magnitude and phase spectra are as expected.
•
Again, we see a performance gain over the previous method and a lower computation time . Due to the power-of-four decimation, fewer recursive iterations are needed which improves the computational speed.
•
The Radix-4 FFT algorithm was also implemented by closely following the theory in order to achieve the most straightforward computation. This can be seen as a modified Radix-2 FFT when looking at the actual code implementation.
70
4.6
Conclusion
In conclusion, this part of the project was useful in demonstrating various computational improvements in FFT algorithms. •
Any method of FFT computation is extremely faster than the mathematical definition implementation of the DFT.
•
These improvements in speed come first from the removal of repeated computation in exponential terms around the unit circle. Secondly, the signals may be decimated in order to further remove repeated calculations.
•
Matlab’s FFT algorithm is many times faster than even our best FFT, showing a possible difference in the way it is handled in the Matlab environment.
71
4.7 • • • • •
•
References Professor SOS, DSP Class notes, University of Texas at San Antonio, Fall 2004. Oppenheim, Schafer, Buck, Discrete-Time Signal Processing, Second Edition, Prentice Hall 1999. Michael J. Andrews, Implementation and Comparison of Radix-2 and Radix-4 FFT Algorithms, Convergent Technology Center, ECE Department, WPI , Worcester, MA 01609-2280 Dr. R. L. Herman, Fourier Analysis of Time Series, UNC Wilmington. L. ROBERT MORRIS,TIME-EFFICIENT RADIX-4 FAST FOURIER TRANSFORM, DEPARTMENT OF SYSTEMS ENGINEERING AND COMPUTING SCIENCE,CARLETON UNIVERSITY, OTTAWA, CANADA MJ Olesen, multivariate complex Fourier transform Queen's University at Kingston, 1995-97
72
4.8
Matlab Code
4.8.1 The code for the main for the radix-2 algorithm % Radix-2 FFT (radix2) % This is the file were the input function is generated % and the Radix-2 function is called. The plotting and time calculation is % done in this main. Adel Hussein clear, clc % Clear variables and screen A = 64; % Period length B = A/2; % Pulse length N = 64*A; % Total duration % Frequency scale for plotting w = (-2*pi:2*pi/N:2*pi-1/N); % Generate one pulse x = [ones(1,B) zeros(1,A-B)]; % Form the repeated pulse chain for count = 0:5 x = [x x]; end % This is to illustrate the correctness using the example in the notes by % professor SOS %x=[4,-3,2,0,-1,-2,3,1]; % This is to calclate the CPU time needed to execute this routine. time = 0; time = cputime; % Initial/start time X = radix2_adelh(x); % Take the Radix-2 FFT time = cputime - time % Execution time % Round the transform for phase X = round(X.*10^6)./(10^6); X = [X X]/N; % Arrange for plotting figure(2) % Plot the magnitude and phase of the inputs FFT subplot(2,1,1) stem(w,abs(X),'.k'), axis([-2*pi 2*pi 0 .55]) title('Radix-2 FFT Magnitude Spectrum of X(\omega)') xlabel('Frequency in Radians \omega'),ylabel('Magnitude') subplot(2,1,2) stem(w,pangle_adelh(X),'.k'), axis([-2*pi 2*pi -2 2])
73
title('Radix-2 FFT Phase Spectrum of X(\omega)') xlabel('Frequency in Radians \omega'),ylabel('Phase (rad)')
4.8.2 The function for the radix-2 algorithm
function X = radix2_adelh(x) % This function recursively finds the % FFT through radix-2 methodology N = length(x); % Find length N of the input data if N == 1 % Test for base case X = x; % Recursive base case else % Apply recursion N2 = N/2; % Transform size WN = cos(2*pi/N)-j*sin(2*pi/N) %exp(-j*2*pi/N); W for computation x1 = x(1:2:end); % Even terms x2 = x(2:2:end); % Odd terms G1 = radix2_adelh(x1); % Recursive call with even terms G2 = radix2_adelh(x2); % Recursive call with odd terms H2 = []; % Initialize array for k = 0:N2-1 % Compute N/2-1 point DFT H2(k+1) = G2(k+1)*WN^k; % Find H2 X(k+1) = G1(k+1)+H2(k+1); % G1+H2 X(k+1+N2) = G1(k+1)-H2(k+1);% G1-H2 end end
4.8.3: The code for the main for the radix-4 algorithm % This is part of project 2. This will calculate FFT using the radix 4 % techniques. % This is the file were the input function is generated % and the Radix-2 function is called. The plotting and time calculation is % done in this main. Adel Hussein % Radix-4 FFT (radix4_adelh) %Initial parameters for the input pulse. clear, clc % Clear variables and screen A = 64; % Period length B = A/2; % Pulse length
74
N = 64*A; % Total duration % Setup the Frequency scale for plotting w = (-2*pi:2*pi/N:2*pi-1/N); % Generate one pulse x = [ones(1,B) zeros(1,A-B)]; % Form the repeated pulse chain for count = 0:5 x = [x x]; end % This is to test the example that was given in the class. Since the example has % only 8 points in the input, it was padded with zeros at the end %x=[4,-3,2,0,-1,-2,3,1,0,0,0,0,0,0,0,0]; % This is to calclate the CPU time needed to execute this routine. N = length(x) % Find length N time = 0; time = cputime; % Initial/start time X = radix4_adelh(x); % Take the Radix-4 FFT time = cputime - time % Execution time
%no_add, no_mult %cmplx1=3*N/8*log2(N); %cmplx2=N*log2(N); %no_add; %no_mult; %fprintf('-----------------\n'); %fprintf('(3N/8)Log(N), NLog(N) , No of Adds , No of Mult\n'); %fprintf('%g , %g , %g , %g\n',cmplx1,cmplx2,no_add,no_mult); %fprintf('-----------------\n'); % Round the transform for phase X = round(X.*10^6)./(10^6); X = [X X]/N; % Arrange for plotting figure(1) % Plot the magnitude and phase subplot(2,1,1) stem(w,abs(X),'.k') axis([-2*pi 2*pi 0 .55]) title(' Radix-4 FFT Magnitude Spectrum of X(\omega)') xlabel('Frequency in Radians \omega');
75
ylabel('Magnitude') subplot(2,1,2) stem(w,pangle_adelh(X),'.k') axis([-2*pi 2*pi -2 2]) title(' Radix-4 FFT Phase Spectrum of X(\omega)') xlabel('Frequency in Radians \omega') ylabel('Phase (rad)') % testing the radix4 function cmpared to the Matlab fft function fs = 32; t = -4:1/fs:4-1/fs; ts = 0.5:1/fs:0.5-1/fs; s1 = (3*cos(2*pi*3*t) + 1*cos(2*pi*5*t))/4 + 5*cos(2*pi*7*t); radix4_fft= radix4_adelh(s1); matlab_fft= fft(s1); nu_samples = 256; freq = -fs/2:fs/nu_samples:fs/2-1/nu_samples; s1_ifft_radix4 = ifft(radix4_fft, length (s1)); s1_error = s1 - s1_ifft_radix4; s1_error_max = max (abs(s1_error)); s1_error_rms = std (s1_error); s1_snr = 10*log10 ( (std(s1))^2/ (std(s1_error))^2); out_s1 = sprintf(' ERROR Measurements, Error = S1 - S1_ifft_radix4 \n'); out_s2 = sprintf(' Error-Max = %g \n ', s1_error_max); out_s3 = sprintf(' Error-RMS = %g \n ', s1_error_rms); out_s4 = sprintf(' SNR = %g \n ', s1_snr); Errror_measures = sprintf('%s %s %s %s',out_s1, out_s2, out_s3, out_s4); display (Errror_measures);
figure (2); subplot(4,1,1); plot (t, s1), title ('Input Signal '); subplot(4,1,2); plot (freq, abs(fftshift(radix4_fft))), title ('Radix FTT Results '); subplot(4,1,3); plot (t, s1_ifft_radix4), title ('IFFT of Radix4 FFT '); subplot(4,1,4); plot (t, s1_error),
76
title ('Error ') xlabel('Time');
fft_error = matlab_fft - radix4_fft; fft_error_max = max (abs(fft_error)); fft_error_rms = std (fft_error); fft_snr = 10*log10 ( (std(matlab_fft))^2/ (std(fft_error))^2); out_ss1 = sprintf(' ERROR Measurements, Error = Matlab-FFT - Radix4-FFT \n'); out_ss2 = sprintf(' Error-Max = %g \n ', fft_error_max); out_ss3 = sprintf(' Error-RMS = %g \n ', fft_error_rms); out_ss4 = sprintf(' SNR = %g \n ', fft_snr); Errror_measures = sprintf('%s %s %s %s',out_ss1, out_ss2, out_ss3, out_ss4); display (Errror_measures);
4.8.4 The function for the radix- algorithm function X = radix4_adelh(x) % This function recursively finds the % FFT through radix-4 methodology % initialize the add and multiply counters.
N = length(x); % Find length N if N == 1 % Test for base case X = x; % Recursive base case else % Apply recursion N4 = N/4; % Transform size W2 = cos(2*pi/N)-j*sin(2*pi/N); %exp(-j*2*pi/N); % Find the 3 WN terms W3 = cos(2*pi*2/N)-j*sin(2*pi*2/N); %exp(-j*2*pi*2/N); % needed to calculate W4 = cos(2*pi*3/N)-j*sin(2*pi*3/N); %exp(-j*2*pi*3/N); % the 4-point DFTs x1 = x(1:4:end); % x(4n) terms x2 = x(2:4:end); % x(4n+1) terms x3 = x(3:4:end); % x(4n+2) terms x4 = x(4:4:end); % x(4n+3) terms G1 = radix4_adelh(x1); % Recursive call using the first group G2 = radix4_adelh(x2); % Recursive call using the second group G3 = radix4_adelh(x3); % Recursive call using the third group G4 = radix4_adelh(x4); % Recursive call using the fourth group H2 = []; % Initialize array for k = 0:N4-1 % Compute N/4-1 point DFT
77
H2(k+1) = G2(k+1)*W2^k; % Find G2 H3(k+1) = G3(k+1)*W3^k; % Find G3 H4(k+1) = G4(k+1)*W4^k; % Find G4 % Apply matrix [X(0,q) [1 1 1 1 [W0*G(0,q) % X(1,q) 1 -j -1 j W1*G(1,q) % X(2,q) = 1 -1 1 -1 W2*G(2,q) % X(3,q)] 1 j -1 -j] W3*G(3,q)] X(k+1) = G1(k+1)+H2(k+1)+H3(k+1)+H4(k+1); X(k+1+N4) = G1(k+1)-j*H2(k+1)-H3(k+1)+j*H4(k+1); X(k+1+2*N4) = G1(k+1)-H2(k+1)+H3(k+1)-H4(k+1); X(k+1+3*N4) = G1(k+1)+j*H2(k+1)-H3(k+1)-j*H4(k+1); end end
4.8.5 The Phase angle Function Code function phase_angle = pangle_adelh(M) % PANGLE: This function correctly computes the 4 % quadrant arctangent of Fourier signal M reM=zeros(1,length(M)); imM=zeros(1,length(M)); reM = real(M); imM = imag(M); phase_angle = zeros(1,length(M)); for h = 1:length(M) if (sign(reM(h))==1 & sign(imM(h))==1) phase_angle(h) = atan(imM(h)/reM(h)); elseif (sign(reM(h))==1 & sign(imM(h))==-1) phase_angle(h) = atan(imM(h)/reM(h)); elseif (sign(reM(h))==-1 & sign(imM(h))==-1) phase_angle(h) = atan(imM(h)/reM(h))-pi; elseif(sign(reM(h))==-1 & sign(imM(h))==1) phase_angle(h) = atan(imM(h)/reM(h))+pi; elseif(sign(reM(h))==0 & sign(imM(h))==0) phase_angle(h) = 0; elseif(sign(reM(h))==0 & sign(imM(h))==1) phase_angle(h) = pi/2; elseif(sign(reM(h))==0 & sign(imM(h))==-1) phase_angle(h) = -pi/2; elseif(sign(reM(h))==1 & sign(imM(h))==0) phase_angle(h) = 0; elseif(sign(reM(h))==-1 & sign(imM(h))==0) phase_angle(h) = pi; end end
78
5 Fast Fourier Transform: Matrix Form (Bassam Jamil Module) Fast Fourier Transform can be computed using matrix representation. The following section explains the theory and implement of matrix representation. We base our discussion mostly on N=8, then we present formulas for general case N. Finally, we present our simulation and analysis results.
5.1 FFT Matrix-Form Theory The Fourier Transform Matrix A is defined as: F[n] = nA* f[n] ……… (1) Let us assume that n=8. The Matrix A can be represented as
(2)
where Wnk = e -2*pi*i*k/n We can make the following observation about the FFT matrix. 1. Wnk is periodic in k and the period is n. Therefore, we do not need to evaluate all the powers between W11 and W(n-1)(n-1). It can be proven that in order to evaluate the FFT matrix, all we need to evaluate Wn1 … W(n)(n-1). However, the challenge is to insert those terms cleverly in the appropriate places in the matrix. 2. The matrix is symmetric Let us shuffle the entries of the vector f[k]: [1, 2, 3, 4, 5, 6, 7, 8] [1, 3, 5, 7, 2, 4, 6, 8] Equation (1) can be rewritten as:
(3)
79
We can observe that: (W81, W83, W85, W87) (W82, W86, W82, W86) (W83, W81, W87, W85) Also, we can observe that:
= W81. (1, W82, W84, W86) = W82. (1, W84, 1 , W84) = W83. (1, W86, W84, W82)
…. (4)
(5)
Equation (5) represents the FFT matrix for 4-point. Combing equations (3) and (4) we obtain:
(6)
Let us examine the lower right corner of 8A (7)
where W81 = -1. Therefore, the lower part of matrix 8A is simply - 8D4. 4A Hence, equation (3) can be expressed as: 8
F(f[1:8]) =
=
=
(8)
Using same argument, we can formulate 4F(f[1:8:2]) and 4F(f[2:8:2]) as:
80
4
F(f[1:8:2]) =
4
F(f[2:8:2]) =
(9)
We can repeat the process for one more iteration: 2
F(f[1:8:4]) =
2
F(f[3:8:4]) =
2
F(f[2:8:4]) =
2
F(f[4:8:4]) =
(10)
where 1F(fk)= fk
5.2 Computation of FFT Matrix-Form The computation of the FFT-Matrix form takes place in layers of matrix computations, as shown in Figure 21, for N=8. The layers are: • Layer-0 reshuffles the input by multiplying input x[n] with the permutation matrix (P). Figure 22 shows the permutation matrix for N=8. In general the permutation matrix can be described as: a. All of the entries are zeros, except the entry P(i, bit-reversal of i) = 1 • The next layers recursively generates the F matrices, where a. Fn = An * F(n/2) b. Figure 23 shows the A’s Matrices.
81
Figure 21 Layers of Matrix Computation
82
Figure 22 Permutation Matrix (N=8) 5.2.1 A’s and D’s Matrices As shown in Figure 23, the A matrix is consist of 4 quadrants. • The left quadrants are the unit matrix of size n/2. • The upper right quadrant is the nDn/2 matrix. This matrix is a diagonal matrix, which all entries are zeros except the diagonal: 1, Wn1, Wn2… Wn(n/2 -1)
83
Figure 23 A's and D's Matrices
5.3 MatLab program Flow Figure 24 shows the Matlab program flow. It consists mainly of the following tasks: • Task1 : Reshuffle the input points by multiplying the input with permutation matrix
84
•
Task2: Generate the output result by looping in two nested loops o The outer loop goes over the layers (see Figure 21) o The inner loop generates the result of the Layer(i) Observe that the complexity of the inner loop is effective O(n). Also, the outer loop is O(log n). Therefore, the complexity of the program is O(n * log n ).
85
86
Figure 24 Program Flow
5.4 Results 5.4.1 4-point Example from Professor Sos Notes This section details the results of the program using the 4-point FFT example in Dr Sos notes, where x = [5, 0, -3, 4] Dr Sos’s note solution is: X [0] = 6 X [1] = 8 + 4*j X [2] = -2 X [3] = 8 – 4*j The FFT-Matrix form program results are (note in shift in Matlab indices.) Sos_4_point_fft = 6.0000 4.0000i
8.0000 + 4.0000i
-2.0000
8.0000 -
Clearly, the program results match the notes solution.
5.4.2 8-point Example from Professor Sos Notes This section details the results of the program using the 8-point FFT example in Dr Sos notes, where x = [ 4, -3, 2, 0, -1, -2, 3, 1] Dr Sos’s note solution is: X [0] = 4 X [1] = 5+j + j * sqrt(2) X [2] = -2 + 6*j X [3] = 5 – j + j*sqrt(2) X [4] = 12 X [5] = 5+j - j * sqrt(2) X [6] = -2 - 6*j X [7] = 5 – j - j*sqrt(2) The FFT-Matrix form program results are (note in shift in Matlab indices.) Sos_8_point_fft = Columns 1 through 5
87
4.0000 0.4142i 12.0000
5.0000 + 2.4142i
-2.0000 + 6.0000i
5.0000 +
Columns 6 through 8 5.0000 - 0.4142i
-2.0000 -6.0000i
5.0000 - 2.4142i
Clearly, the program results match the notes solution.
5.4.3 Compare Matrix Form FFT with Matlab Functions In this section we compare Matrix Form FFT with matlab functions. We will try to address the following issues: • How accurate is the reconstructed signal compared to original signal, • How accurate is the results of Matrix-Form FFT to matlab FFT • How fast is the Matrix-Form FFT to matlab FFT I. Input Signal The input signal used in the calculations is: x(t) = 3*cos(2*pi*3*t) + 1*cos(2*pi*5*t) + 5*cos(2*pi*7*t); (See Figure 25, the first graph from top) The above signal will be sampled at FS=16, between t=-4 and t=-4-1/FS. Therefore , the discrete x[n] is sampled from signal x(t) as follows: x[1] = x(-4.0000), x[2] = x(-3.9375), x[3] = x(-3.8750) … x[127] = x(3.8750), x[128] = x(3.9375); II. Matrix-Form FFT Next, we calculate Matrix-Form-FFT. Figure 25 ( second graph from top) shows a plot for the frequency spectrum of the matrix-form FFT. Notice that the frequency components are showing up at • Frequency=3 • Frequency=5 • Frequency=7 This matches the x(t) signal frequencies. III. Next, we calculate the IFFT of the Matrix-Form FFT We reconstruct the signal using IFFT function, see reconstructed signal in Figure 25 (third graph from top). We them calculate the ERROR ERROR = (original sampled signal) - (IFFT of the Matrix-Form FFT) See error signal in Figure 25 (the bottom graph)
88
Error results are:
ERROR Type Maximum Error ERROR RMS SNR
Value 1.60119e-015 7.29786e-016 315.201
IV. Finally, we calculate FFT using Matlab FFT function We then calculate ERROR ERROR = (FFT results from Matlab) – ( FFT results from Matrix-Form Algorithm) Error results are: ERROR Type Value Maximum Error 5.16749e-014 ERROR RMS 6.78428e-015 V. Finally, we compare cputime Next, we will measure the CPU Time. For that we will increase the range of time to: • fs=16 • t = -64:1/fs:64-1/fs;
FFT Type CPUTIME (Second) Matlab 0.010000 Matrix Form 2.033000 Note: To calculate cputime, we need to disable the reporting of number of operations in the FFT Matrix Form
Conclusions: • The matrix-form computations are very accurate. The SNR is excellent. • FFT Matrix Form is slower. There are many areas of improvement. Example o The multiplication with permutation matrix is not needed. Instead, simple bit-reversal does the job.
89
90
Figure 25 FFT: Matrix Form vs. Matlab 5.4.4 Number of operations In this section we calculate number of operations in the FFT matrix form algorithm. First let us start with the following definitions: Terminology Generic Operation
Definition Mathematical operation between 2 non zero complex numbers Addition of Real part or Imaginary part of complex numbers Multiplication of Complex number components (Real or Imaginary)
ADD operation MULT operation
Example: For the operation (5+9j) * ( 6+8j) • It is considered 1 Generic operation • It has 4 MULTs • It has 2 Additions The rest of the section explores number of operations in 8-point FFT Matrix-Form with different inputs.
Case of x[n] [0,0,0,0, 0,0,0,0] [1,1,1,1, 0,0,0,0] [0,0,0,0, 1,1,1,1] [1,1,1,1, 1,1,1,1] [1,0,1,0, 1,0,1,0] [0,1,0, 1,0,1,0,1] Non (1 or 0)
ADD Op 0 8 8 0 0 0 8
MUL Op 0 9 19 3 1 2 32
GEN Op 0 36 39 28 14 14 48
Conclusion: Number of operations is sensitive to the type of input data. Especially for input data that is rich in 1’s and 0’s. This could be exploited further for performance optimizations.
5.5 Matlab Code The program consists of the following codes: 1. bmohd_fft_matrix_top.m: this is the top-level matlab program. It is used to generate the above results 2. bmohd_fft_matrix_form: this is the matlab code for FFT Matrix-Form
91
3. bmohd_A_kn.m: it contains a function to calculate and A’s D’s matrices 4. bmohd_calc_operations: it calculate the ADD, MULT and Generic operations in complex matrix multiplication.
5.5.1 % % %
Bmohd_fft_matrix_top.m
Author: Bassam Mohd Description: The following tests the FFT Matrix-Form Program
% % Case #1: Calcuate the 4-point FFT example of Dr. Sos Notes ADD=0; MUL=0; GEN=0; Sos_4_point_f = [5, 0, -3, 4] ; Sos_4_point_fft = bmohd_fft_matrix_form (Sos_4_point_f); display (Sos_4_point_fft); % % Case #2: Calcuate the 8-point FFT example of Dr. Sos Notes Sos_8_point_f = [4, -3, 2, 0, -1, -2, 3, 1] ; Sos_8_point_fft = bmohd_fft_matrix_form (Sos_8_point_f); display (Sos_8_point_fft); % Case #3: Compare Matrix-Form FFT resutls with Matlab FFT fs = 16; t = -4:1/fs:4-1/fs; t_small = 0.5:1/fs:0.5-1/fs; x_original = 3*cos(2*pi*3*t) + 1*cos(2*pi*5*t) + 5*cos(2*pi*7*t); matrix_form_fft= bmohd_fft_matrix_form(x_original); matlab_form_fft= fft(x_original); fft_points = 128; freq = -fs/2:fs/fft_points:fs/2-1/fft_points; x_ifft_of_matrix_fft= ifft(matrix_form_fft, length (x_original)); x_error x_error_max x_error_rms x_snr
= x_original - x_ifft_of_matrix_fft; = max (abs(x_error)); = std (x_error); = 10*log10 ( (std(x_original))^2/ (std(x_error))^2);
out_s1 = sprintf(' ERROR Measurements, Error = Original-Signal IFFT-of-Matrix-Form-FFT \n'); out_s2 = sprintf(' Error-Max = %g \n ', x_error_max); out_s3 = sprintf(' Error-RMS = %g \n ', x_error_rms); out_s4 = sprintf(' SNR = %g \n ', x_snr); Errror_measures = sprintf('%s %s %s %s',out_s1, out_s2, out_s3, out_s4); display (Errror_measures); figure (1);
92
subplot(4,1,1); plot (t, x_original), title ('x '); subplot(4,1,2); stem (freq, abs(fftshift(matrix_form_fft))), title ('Matrix-FFT '); subplot(4,1,3); plot (t, x_ifft_of_matrix_fft), title ('IFFT'); subplot(4,1,4); plot (t, x_error), title ('Error '); fft_error = matlab_form_fft - matrix_form_fft; fft_error_max = max (abs(fft_error)); fft_error_rms = std (fft_error); fft_snr = 10*log10 ( (std(matlab_form_fft))^2/ (std(fft_error))^2); out_ss1 = sprintf(' ERROR Measurements, Error = Matlab-FFT - MatrixForm-FFT \n'); out_ss2 = sprintf(' Error-Max = %g \n ', fft_error_max); out_ss3 = sprintf(' Error-RMS = %g \n ', fft_error_rms); out_ss4 = sprintf(' SNR = %g \n ', fft_snr); Errror_measures = sprintf('%s %s %s %s',out_ss1, out_ss2, out_ss3, out_ss4); display (Errror_measures); % Case3.2: Calculate CPU Time fs = 16; t = -64:1/fs:64-1/fs; t_small = 0.5:1/fs:0.5-1/fs; x_original = 3*cos(2*pi*3*t) + 1*cos(2*pi*5*t) + 5*cos(2*pi*7*t); cputime1 = cputime; matrix_form_fft= bmohd_fft_matrix_form(x_original); cputime2 = cputime; matlab_form_fft= fft(x_original); cputime3 = cputime; cputime_s = sprintf(' CPUTIME: Matrix-Form FFT=%f Matlab FFT=%f \n',cputime2-cputime1, cputime3-cputime2); display (cputime_s); % Case #5: Number of Operations % All zeros case4_f = [0, 0, 0, 0, 0, 0, 0, 0] ; [case4_fft ADD, MUL, GEN]= bmohd_fft_matrix_form (case4_f); op_cnt = sprintf(' Case 4.1 : Number of operation ADD=%d MUL =%d GEN=%d', ADD, MUL, GEN); display(op_cnt); % Half Ones case4_f = [1, 1, 1, 1, 0, 0, 0, 0] ; [case4_fft, ADD, MUL, GEN]= bmohd_fft_matrix_form (case4_f); op_cnt = sprintf(' Case 4.2 : Number of operation ADD=%d MUL =%d GEN=%d', ADD, MUL, GEN); display(op_cnt);
93
% HALF Ones case4_f = [0, 0, 0, 0, 1, 1, 1, 1] ; [case4_fft, ADD, MUL, GEN]= bmohd_fft_matrix_form (case4_f); op_cnt = sprintf(' Case 4.3 : Number of operation ADD=%d MUL GEN=%d', ADD, MUL, GEN); display(op_cnt); % ALL Ones case4_f = [1, 1, 1, 1, 1, 1, 1, 1] ; [case4_fft, ADD, MUL, GEN]= bmohd_fft_matrix_form (case4_f); op_cnt = sprintf(' Case 4.4 : Number of operation ADD=%d MUL GEN=%d', ADD, MUL, GEN); display(op_cnt); % Ones/Zeros case4_f = [1, 0, 1, 0, 1, 0, 1, 0] ; [case4_fft, ADD, MUL, GEN]= bmohd_fft_matrix_form (case4_f); op_cnt = sprintf(' Case 4.5 : Number of operation ADD=%d MUL GEN=%d', ADD, MUL, GEN); display(op_cnt); % Ones/Zeros case4_f = [0, 1, 0, 1, 0, 1, 0, 1] ; [case4_fft, ADD, MUL, GEN]= bmohd_fft_matrix_form (case4_f); op_cnt = sprintf(' Case 4.6 : Number of operation ADD=%d MUL GEN=%d', ADD, MUL, GEN); display(op_cnt); % Real > 1 case4_f = [2, 3, 4, 5, 6, 7, 8, 9] ; [case4_fft, ADD, MUL, GEN]= bmohd_fft_matrix_form (case4_f); op_cnt = sprintf(' Case 4.7 : Number of operation ADD=%d MUL GEN=%d', ADD, MUL, GEN); display(op_cnt);
5.5.2 % % %
=%d
=%d
=%d
=%d
=%d
bmohd_fft_matrix_form
Author: Bassam Mohd Description: The following program implement FFT Matrix-Form
% Function Call definition function [ Matrix_FFT, Total_ADD, Total_MUL, Total_GEN ] = bmohd_fft_matrix_form (f); Total_ADD = 0; Total_MUL = 0; Total_GEN = 0; ADD =0; MUL =0; GEN =0; % Calculate the length of input signal N =length (f); % Construct the Permutation Matrix bits = log2 (N); PERM = zeros (N, N); for n=1: N
94
rev_n = 0; for bit_i=1: bits rev_n = rev_n + 2^(bits - bit_i )* bitget(n-1, bit_i); end PERM(rev_n+1,n ) = 1; end % Shuffle the input by multiplying input signal with permutation matrix f_shuffled = f * PERM; %Construct the FFT Matrices F_prev = f_shuffled; % Outer Loop for Layer for layer = 1 : log2(N) N_point = 2^layer; F_nxt = []; %Inner Loop for Sections of a Layer for section = 1 : N/N_point F_A = zeros(N_point, 1); %Construct F sub-marices for q =1 : N_point F_A(q) = F_prev ((section-1)*N_point +q); end %q A_nxt = bmohd_A_kn( N_point/2, N_point); F_tmp = A_nxt * F_A; %Calculation number of operations [GEN, ADD, MUL] = bmohd_calc_operations (A_nxt, F_A); Total_ADD = Total_ADD + ADD; Total_MUL = Total_MUL + MUL; Total_GEN = Total_GEN + GEN; %Construct full 'F' Matrix from Sub-matrices for q =1 : N_point F_nxt((section-1)*N_point +q) = F_tmp(q); end %q end %r F_prev = []; F_prev = F_nxt; end %s % Final Result Matrix_FFT = F_prev; % This Part of the program is used for Debug purposes. I will flag an error % when the program results mismatch Matlab's. MATLAB_FFT = fft(f, N); err_cnt = 0; for r=1: length(MATLAB_FFT) if ( abs (abs(MATLAB_FFT(r)) - abs(Matrix_FFT(r))) > 0.000001) err_cnt = err_cnt +1; end end if ( err_cnt>0) Error_Message = '*** Error***: Matrix Form has greater than 0.000001 mismatch with Matlab.'
95
display (Error_Message); end
5.5.3 % % %
bmohd_A_kn.m
Author: Bassam Mohd Description: This program calcuates nxn A matrix
function [ AA] = bmohd_A_kn (k, N); %% First : Calculate D % Initialize the D matrix D = eye (N/2); % Plug in the W's in the diagonal of D matrix for k=2:N/2 %D(i,i) = ww( i-1, N); D(k,k) = exp (-2*pi*i*(k-1)/N); end %% Second: Construct the A matrix AA = [ eye(N/2) , D; eye(N/2) , -D ];
5.5.4 % % %
bmohd_calc_operations
Author: Bassam Mohd Description: This program calcuates number of operations
function [ Total_GEN, Total_ADD, Total_MUL] = bmohd_calc_operations (A, B); [A_1, A_2] = size(A); [B_1, B_2] = size(B); Total_ADD = 0; Total_MUL = 0; Total_GEN = 0; %%% Loop through matrices and calc operations for a1 =1:A_1 for b1 =1:B_2 for c1 =1:A_2 ADD_op = 2; MUL_op = 4; GEN_op = 1; A(a1,c1); B(c1,b1); a = real(A(a1,c1)); b = imag(A(a1,c1)); c = real(B(c1,b1)); d = imag(B(c1,b1)); if ( a*c == 0 | a == 1 | c == 1)
96
MUL_op = end if ( b*d == MUL_op = end if ( a*c == ADD_op = end
MUL_op -1; 0 | b == 1 | d == 1) MUL_op -1; 0 | b*d==0) ADD_op -1;
if ( a*d == 0 | a == 1 | d == 1) MUL_op = MUL_op -1; end if ( b*c == 0 | b == 1 | c == 1) MUL_op = MUL_op -1; end if ( a*d == 0 | b*c==0) ADD_op = ADD_op -1; end if ( (a==0 & b==0) | (c==0 & d==0) ) GEN_op =0; end Total_ADD = Total_ADD + ADD_op; Total_MUL = Total_MUL + MUL_op; Total_GEN = Total_GEN + GEN_op; end end end
97
