diff --git a/src/main/java/com/maths/CircularConvolutionFFT.java b/src/main/java/com/maths/CircularConvolutionFFT.java new file mode 100644 index 00000000..9182de9e --- /dev/null +++ b/src/main/java/com/maths/CircularConvolutionFFT.java @@ -0,0 +1,58 @@ +package com.maths; + +import java.util.ArrayList; + +/** + * Class for circular convolution of two discrete signals using the convolution theorem. + * + * @author Ioannis Karavitsis + * @version 1.0 + * */ +public class CircularConvolutionFFT +{ + /** + * This method pads the signal with zeros until it reaches the new size. + * + * @param x The signal to be padded. + * @param newSize The new size of the signal. + * */ + private static void padding(ArrayList x, int newSize) + { + if(x.size() < newSize) + { + int diff = newSize - x.size(); + for(int i = 0; i < diff; i++) + x.add(new FFT.Complex()); + } + } + + /** + * Discrete circular convolution function. It uses the convolution theorem for discrete signals: convolved = IDFT(DFT(a)*DFT(b)). + * Then we use the FFT algorithm for faster calculations of the two DFTs and the final IDFT. + * + * More info: + * https://en.wikipedia.org/wiki/Convolution_theorem + * + * @param a The first signal. + * @param b The other signal. + * @return The convolved signal. + * */ + public static ArrayList fftCircularConvolution(ArrayList a, ArrayList b) + { + int convolvedSize = Math.max(a.size(), b.size()); //The two signals must have the same size equal to the bigger one + padding(a, convolvedSize); //Zero padding the smaller signal + padding(b, convolvedSize); + + /* Find the FFTs of both signal. Here we use the Bluestein algorithm because we want the FFT to have the same length with the signal and not bigger */ + FFTBluestein.fftBluestein(a, false); + FFTBluestein.fftBluestein(b, false); + ArrayList convolved = new ArrayList<>(); + + for(int i = 0; i < a.size(); i++) + convolved.add(a.get(i).multiply(b.get(i))); //FFT(a)*FFT(b) + + FFTBluestein.fftBluestein(convolved, true); //IFFT + + return convolved; + } +} diff --git a/src/main/java/com/maths/Convolution.java b/src/main/java/com/maths/Convolution.java new file mode 100644 index 00000000..575c1214 --- /dev/null +++ b/src/main/java/com/maths/Convolution.java @@ -0,0 +1,47 @@ +package com.maths; + +/** + * Class for linear convolution of two discrete signals + * + * @author Ioannis Karavitsis + * @version 1.0 + * */ +public class Convolution +{ + /** + * Discrete linear convolution function. Both input signals and the output signal must start from 0. + * If you have a signal that has values before 0 then shift it to start from 0. + * + * @param A The first discrete signal + * @param B The second discrete signal + * @return The convolved signal + * */ + public static double[] convolution(double[] A, double[] B) + { + double[] convolved = new double[A.length + B.length - 1]; + + /* + The discrete convolution of two signals A and B is defined as: + + A.length + C[i] = Σ (A[k]*B[i-k]) + k=0 + + It's obvious that: 0 <= k <= A.length , 0 <= i <= A.length + B.length - 2 and 0 <= i-k <= B.length - 1 + From the last inequality we get that: i - B.length + 1 <= k <= i and thus we get the conditions below. + */ + for(int i = 0; i < convolved.length; i++) + { + convolved[i] = 0; + int k = Math.max(i - B.length + 1, 0); + + while(k < i + 1 && k < A.length) + { + convolved[i] += A[k] * B[i - k]; + k++; + } + } + + return convolved; + } +} \ No newline at end of file diff --git a/src/main/java/com/maths/ConvolutionFFT.java b/src/main/java/com/maths/ConvolutionFFT.java new file mode 100644 index 00000000..3897be6c --- /dev/null +++ b/src/main/java/com/maths/ConvolutionFFT.java @@ -0,0 +1,61 @@ +package com.maths; + +import java.util.ArrayList; + +/** + * Class for linear convolution of two discrete signals using the convolution theorem. + * + * @author Ioannis Karavitsis + * @version 1.0 + * */ +public class ConvolutionFFT +{ + /** + * This method pads the signal with zeros until it reaches the new size. + * + * @param x The signal to be padded. + * @param newSize The new size of the signal. + * */ + private static void padding(ArrayList x, int newSize) + { + if(x.size() < newSize) + { + int diff = newSize - x.size(); + for(int i = 0; i < diff; i++) + x.add(new FFT.Complex()); + } + } + + /** + * Discrete linear convolution function. It uses the convolution theorem for discrete signals convolved: = IDFT(DFT(a)*DFT(b)). + * This is true for circular convolution. In order to get the linear convolution of the two signals we first pad the two signals to have the same size equal to the convolved signal (a.size() + b.size() - 1). + * Then we use the FFT algorithm for faster calculations of the two DFTs and the final IDFT. + * + * More info: + * https://en.wikipedia.org/wiki/Convolution_theorem + * https://ccrma.stanford.edu/~jos/ReviewFourier/FFT_Convolution.html + * + * @param a The first signal. + * @param b The other signal. + * @return The convolved signal. + * */ + public static ArrayList convolutionFFT(ArrayList a, ArrayList b) + { + int convolvedSize = a.size() + b.size() - 1; //The size of the convolved signal + padding(a, convolvedSize); //Zero padding both signals + padding(b, convolvedSize); + + /* Find the FFTs of both signals (Note that the size of the FFTs will be bigger than the convolvedSize because of the extra zero padding in FFT algorithm) */ + FFT.fft(a, false); + FFT.fft(b, false); + ArrayList convolved = new ArrayList<>(); + + for(int i = 0; i < a.size(); i++) + convolved.add(a.get(i).multiply(b.get(i))); //FFT(a)*FFT(b) + + FFT.fft(convolved, true); //IFFT + convolved.subList(convolvedSize, convolved.size()).clear(); //Remove the remaining zeros after the convolvedSize. These extra zeros came from paddingPowerOfTwo() method inside the fft() method. + + return convolved; + } +} diff --git a/src/main/java/com/maths/FFT.java b/src/main/java/com/maths/FFT.java new file mode 100644 index 00000000..6a8f1217 --- /dev/null +++ b/src/main/java/com/maths/FFT.java @@ -0,0 +1,281 @@ +package com.maths; + +import java.util.ArrayList; +import java.util.Collections; + +/** + * Class for calculating the Fast Fourier Transform (FFT) of a discrete signal using the Cooley-Tukey algorithm. + * + * @author Ioannis Karavitsis + * @version 1.0 + * */ +public class FFT +{ + /** + * This class represents a complex number and has methods for basic operations. + * + * More info: + * https://introcs.cs.princeton.edu/java/32class/Complex.java.html + * */ + static class Complex + { + private double real, img; + + /** + * Default Constructor. + * Creates the complex number 0. + * */ + public Complex() + { + real = 0; + img = 0; + } + + /** + * Constructor. Creates a complex number. + * + * @param r The real part of the number. + * @param i The imaginary part of the number. + * */ + public Complex(double r, double i) + { + real = r; + img = i; + } + + /** + * Returns the real part of the complex number. + * + * @return The real part of the complex number. + * */ + public double getReal() + { + return real; + } + + /** + * Returns the imaginary part of the complex number. + * + * @return The imaginary part of the complex number. + * */ + public double getImaginary() + { + return img; + } + + /** + * Adds this complex number to another. + * + * @param z The number to be added. + * @return The sum. + * */ + public Complex add(Complex z) + { + Complex temp = new Complex(); + temp.real = this.real + z.real; + temp.img = this.img + z.img; + return temp; + } + + /** + * Subtracts a number from this complex number. + * + * @param z The number to be subtracted. + * @return The difference. + * */ + public Complex subtract(Complex z) + { + Complex temp = new Complex(); + temp.real = this.real - z.real; + temp.img = this.img - z.img; + return temp; + } + + /** + * Multiplies this complex number by another. + * + * @param z The number to be multiplied. + * @return The product. + * */ + public Complex multiply(Complex z) + { + Complex temp = new Complex(); + temp.real = this.real*z.real - this.img*z.img; + temp.img = this.real*z.img + this.img*z.real; + return temp; + } + + /** + * Multiplies this complex number by a scalar. + * + * @param n The real number to be multiplied. + * @return The product. + * */ + public Complex multiply(double n) + { + Complex temp = new Complex(); + temp.real = this.real * n; + temp.img = this.img * n; + return temp; + } + + /** + * Finds the conjugate of this complex number. + * + * @return The conjugate. + * */ + public Complex conjugate() + { + Complex temp = new Complex(); + temp.real = this.real; + temp.img = -this.img; + return temp; + } + + /** + * Finds the magnitude of the complex number. + * + * @return The magnitude. + * */ + public double abs() + { + return Math.hypot(this.real, this.img); + } + + /** + * Divides this complex number by another. + * + * @param z The divisor. + * @return The quotient. + * */ + public Complex divide(Complex z) + { + Complex temp = new Complex(); + temp.real = (this.real*z.real + this.img*z.img) / (z.abs()*z.abs()); + temp.img = (this.img*z.real - this.real*z.img) / (z.abs()*z.abs()); + return temp; + } + + /** + * Divides this complex number by a scalar. + * + * @param n The divisor which is a real number. + * @return The quotient. + * */ + public Complex divide(double n) + { + Complex temp = new Complex(); + temp.real = this.real / n; + temp.img = this.img / n; + return temp; + } + } + + /** + * Iterative In-Place Radix-2 Cooley-Tukey Fast Fourier Transform Algorithm with Bit-Reversal. + * The size of the input signal must be a power of 2. If it isn't then it is padded with zeros and the output FFT will be bigger than the input signal. + * + * More info: + * https://www.algorithm-archive.org/contents/cooley_tukey/cooley_tukey.html + * https://www.geeksforgeeks.org/iterative-fast-fourier-transformation-polynomial-multiplication/ + * https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm + * https://cp-algorithms.com/algebra/fft.html + * + * @param x The discrete signal which is then converted to the FFT or the IFFT of signal x. + * @param inverse True if you want to find the inverse FFT. + * */ + public static void fft(ArrayList x, boolean inverse) + { + /* Pad the signal with zeros if necessary */ + paddingPowerOfTwo(x); + int N = x.size(); + + /* Find the log2(N) */ + int log2N = 0; + while((1 << log2N) < N) + log2N++; + + /* Swap the values of the signal with bit-reversal method */ + int reverse; + for(int i = 0; i < N; i++) + { + reverse = reverseBits(i, log2N); + if(i < reverse) + Collections.swap(x, i, reverse); + } + + int direction = inverse ? -1 : 1; + + /* Main loop of the algorithm */ + for(int len = 2; len <= N; len *= 2) + { + double angle = -2 * Math.PI / len * direction; + Complex wlen = new Complex(Math.cos(angle), Math.sin(angle)); + for(int i = 0; i < N; i += len) + { + Complex w = new Complex(1, 0); + for(int j = 0; j < len / 2; j++) + { + Complex u = x.get(i + j); + Complex v = w.multiply(x.get(i + j + len/2)); + x.set(i + j, u.add(v)); + x.set(i + j + len/2, u.subtract(v)); + w = w.multiply(wlen); + } + } + } + + /* Divide by N if we want the inverse FFT */ + if(inverse) + { + for (int i = 0; i < x.size(); i++) + { + Complex z = x.get(i); + x.set(i, z.divide(N)); + } + } + } + + /** + * This function reverses the bits of a number. + * It is used in Cooley-Tukey FFT algorithm. + * + * E.g. + * num = 13 = 00001101 in binary + * log2N = 8 + * Then reversed = 176 = 10110000 in binary + * + * More info: + * https://cp-algorithms.com/algebra/fft.html + * https://www.geeksforgeeks.org/write-an-efficient-c-program-to-reverse-bits-of-a-number/ + * + * @param num The integer you want to reverse its bits. + * @param log2N The number of bits you want to reverse. + * @return The reversed number + * */ + private static int reverseBits(int num, int log2N) + { + int reversed = 0; + for(int i = 0; i < log2N; i++) + { + if((num & (1 << i)) != 0) + reversed |= 1 << (log2N - 1 - i); + } + return reversed; + } + + /** + * This method pads an ArrayList with zeros in order to have a size equal to the next power of two of the previous size. + * + * @param x The ArrayList to be padded. + * */ + private static void paddingPowerOfTwo(ArrayList x) + { + int n = 1; + int oldSize = x.size(); + while(n < oldSize) + n *= 2; + for(int i = 0; i < n - oldSize; i++) + x.add(new Complex()); + } +} diff --git a/src/main/java/com/maths/FFTBluestein.java b/src/main/java/com/maths/FFTBluestein.java new file mode 100644 index 00000000..644295cb --- /dev/null +++ b/src/main/java/com/maths/FFTBluestein.java @@ -0,0 +1,69 @@ +package com.maths; + +import java.util.ArrayList; + +/** + * Class for calculating the Fast Fourier Transform (FFT) of a discrete signal using the Bluestein's algorithm. + * + * @author Ioannis Karavitsis + * @version 1.0 + * */ +public class FFTBluestein +{ + /** + * Bluestein's FFT Algorithm. + * + * More info: + * https://en.wikipedia.org/wiki/Chirp_Z-transform#Bluestein.27s_algorithm + * http://tka4.org/materials/lib/Articles-Books/Numerical%20Algorithms/Hartley_Trasform/Bluestein%27s%20FFT%20algorithm%20-%20Wikipedia,%20the%20free%20encyclopedia.htm + * + * @param x The discrete signal which is then converted to the FFT or the IFFT of signal x. + * @param inverse True if you want to find the inverse FFT. + * */ + public static void fftBluestein(ArrayList x, boolean inverse) + { + int N = x.size(); + int bnSize = 2*N - 1; + int direction = inverse ? -1 : 1; + ArrayList an = new ArrayList<>(); + ArrayList bn = new ArrayList<>(); + + /* Initialization of the b(n) sequence (see Wikipedia's article above for the symbols used)*/ + for(int i = 0; i < bnSize; i++) + bn.add(new FFT.Complex()); + + for(int i = 0; i < N; i++) + { + double angle = (i - N + 1) * (i - N + 1) * Math.PI / N * direction; + bn.set(i, new FFT.Complex(Math.cos(angle), Math.sin(angle))); + bn.set(bnSize - i - 1, new FFT.Complex(Math.cos(angle), Math.sin(angle))); + } + + /* Initialization of the a(n) sequence */ + for(int i = 0; i < N; i++) + { + double angle = -i * i * Math.PI / N * direction; + an.add(x.get(i).multiply(new FFT.Complex(Math.cos(angle), Math.sin(angle)))); + } + + ArrayList convolution = ConvolutionFFT.convolutionFFT(an, bn); + + /* The final multiplication of the convolution with the b*(k) factor */ + for(int i = 0; i < N; i++) + { + double angle = -1 * i * i * Math.PI / N * direction; + FFT.Complex bk = new FFT.Complex(Math.cos(angle), Math.sin(angle)); + x.set(i, bk.multiply(convolution.get(i + N - 1))); + } + + /* Divide by N if we want the inverse FFT */ + if(inverse) + { + for (int i = 0; i < N; i++) + { + FFT.Complex z = x.get(i); + x.set(i, z.divide(N)); + } + } + } +} diff --git a/src/test/java/com/maths/CircularConvolutionFFTTest.java b/src/test/java/com/maths/CircularConvolutionFFTTest.java new file mode 100644 index 00000000..c18922e5 --- /dev/null +++ b/src/test/java/com/maths/CircularConvolutionFFTTest.java @@ -0,0 +1,68 @@ +package com.maths; + +import static org.junit.jupiter.api.Assertions.*; +import org.junit.jupiter.api.Test; + +import java.math.BigDecimal; +import java.math.RoundingMode; +import java.util.ArrayList; + +/** + * Class for testing the circular convolution of two signals. + * You can also use Matlab for bigger signals and compare the results by modifying the code below. + * + * @author Ioannis Karavitsis + * @version 1.0 + * */ +class CircularConvolutionFFTTest +{ + /** + * Function to round a number with double precision to n decimal places. + * + * More info: + * https://www.baeldung.com/java-round-decimal-number + * + * @param num The number to be rounded. + * @param places The number of decimal places up to which the number will be rounded. + * @return The rounded number. + * */ + private double round(double num, int places) + { + if (places < 0) throw new IllegalArgumentException(); + + BigDecimal bd = new BigDecimal(Double.toString(num)); + bd = bd.setScale(places, RoundingMode.HALF_UP); + return bd.doubleValue(); + } + + @Test + public void testCircularConvolutionFFT() + { + ArrayList x = new ArrayList<>(); + for (int i = 1; i <= 10; i++) //x = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 } + x.add(new FFT.Complex(i, 0)); + + ArrayList y = new ArrayList<>(); + for (int i = 0; i < 5; i++) //y = { 0.2, 0.2, 0.2, 0.2, 0.2 } + y.add(new FFT.Complex(0.2, 0)); + + /* Matlab results */ + double[] res = { 7, 6, 5, 4, 3, 4, 5, 6, 7, 8 }; + + /* Find the circular convolution of the two signals */ + ArrayList convolved = CircularConvolutionFFT.fftCircularConvolution(x, y); + + /* Compare the results with those from Matlab (command: cconv(x, y, 10) for other signals you must use the larger size of the two signals, here is 10). + I rounded them to 1 decimal place, due to the different calculations and precision that Matlab and this Java program use. */ + for(int i = 0; i < convolved.size(); i++) + { + assertEquals(res[i], round(convolved.get(i).getReal(), 1)); + assertEquals(0, round(convolved.get(i).getImaginary(), 1)); + } + + /* Print the results (Optional) */ + /*for (FFT.Complex complex : convolved) + System.out.println(complex.getReal() + " " + complex.getImaginary());*/ + } + +} \ No newline at end of file diff --git a/src/test/java/com/maths/ConvolutionFFTTest.java b/src/test/java/com/maths/ConvolutionFFTTest.java new file mode 100644 index 00000000..8a1f4072 --- /dev/null +++ b/src/test/java/com/maths/ConvolutionFFTTest.java @@ -0,0 +1,69 @@ +package com.maths; + +import static org.junit.jupiter.api.Assertions.*; +import org.junit.jupiter.api.Test; +import java.math.BigDecimal; +import java.math.RoundingMode; +import java.util.ArrayList; + +/** + * Class for testing the Convolution of two signals using the Cooley-Tukey FFT algorithm. + * You can also use Matlab for bigger signals and compare the results by modifying the code below. + * + * @author Ioannis Karavitsis + * @version 1.0 + * */ +class ConvolutionFFTTest +{ + /** + * Function to round a number with double precision to n decimal places. + * + * More info: + * https://www.baeldung.com/java-round-decimal-number + * + * @param num The number to be rounded. + * @param places The number of decimal places up to which the number will be rounded. + * @return The rounded number. + * */ + private double round(double num, int places) + { + if (places < 0) throw new IllegalArgumentException(); + + BigDecimal bd = new BigDecimal(Double.toString(num)); + bd = bd.setScale(places, RoundingMode.HALF_UP); + return bd.doubleValue(); + } + + /** + * Test for convolutionFFT() function. + * */ + @Test + public void testConvolutionFFT() + { + ArrayList x = new ArrayList<>(); + for (int i = 1; i <= 10; i++) //x = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 } + x.add(new FFT.Complex(i, 0)); + + ArrayList y = new ArrayList<>(); + for (int i = 0; i < 5; i++) //y = { 0.2, 0.2, 0.2, 0.2, 0.2 } + y.add(new FFT.Complex(0.2, 0)); + + /* Matlab results */ + double[] real = { 0.2, 0.6, 1.2, 2, 3, 4, 5, 6, 7.000000000000001, 8, 6.800000000000001, 5.4, 3.8, 2 }; + + /* Find the convolution of the two signals */ + ArrayList convolved = ConvolutionFFT.convolutionFFT(x, y); + + /* Compare the results with those from Matlab (command: conv(x, y)). + I rounded them to 1 decimal place, due to the different calculations and precision that Matlab and this Java program use. */ + for(int i = 0; i < convolved.size(); i++) + { + assertEquals(round(real[i], 1), round(convolved.get(i).getReal(), 1)); + assertEquals(0, round(convolved.get(i).getImaginary(),1)); + } + + /* Print the results (Optional) */ + /*for (FFT.Complex complex : convolved) + System.out.println(complex.getReal() + " " + complex.getImaginary());*/ + } +} \ No newline at end of file diff --git a/src/test/java/com/maths/ConvolutionTest.java b/src/test/java/com/maths/ConvolutionTest.java new file mode 100644 index 00000000..0ce6aec7 --- /dev/null +++ b/src/test/java/com/maths/ConvolutionTest.java @@ -0,0 +1,56 @@ +package com.maths; + +import static org.junit.jupiter.api.Assertions.*; +import org.junit.jupiter.api.Test; +import java.math.BigDecimal; +import java.math.RoundingMode; + +/** + * Class for testing the convolution of two discrete signals. + * You can also use Matlab for bigger signals and compare the results by modifying the code below. + * + * @author Ioannis Karavitsis + * @version 1.0 + * */ +public class ConvolutionTest +{ + /** + * Function to round a number with double precision to n decimal places. + * + * More info: + * https://www.baeldung.com/java-round-decimal-number + * + * @param num The number to be rounded. + * @param places The number of decimal places up to which the number will be rounded. + * @return The rounded number. + * */ + private double round(double num, int places) + { + if (places < 0) throw new IllegalArgumentException(); + + BigDecimal bd = new BigDecimal(Double.toString(num)); + bd = bd.setScale(places, RoundingMode.HALF_UP); + return bd.doubleValue(); + } + + /** + * Test for convolution() function. + * */ + @Test + public void testConvolution() + { + double[] x = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 }; + double[] y = { 0.2, 0.2, 0.2, 0.2, 0.2 }; + + /* Matlab results */ + double[] res = { 0.2, 0.6, 1.2, 2, 3, 4, 5, 6, 7.000000000000001, 8, 6.800000000000001, 5.4, 3.8, 2 }; + + /* Find the convolution f the two signals */ + double[] convolved = Convolution.convolution(x, y); + + /* Compare the results with those from Matlab (command: conv(x, y)). + I rounded them to 1 decimal place, due to the different calculations and precision that Matlab and this Java program use. */ + for(int i = 0; i < convolved.length; i++) + assertEquals(round(res[i], 1), round(convolved[i], 1)); + } +} \ No newline at end of file diff --git a/src/test/java/com/maths/FFTBluesteinTest.java b/src/test/java/com/maths/FFTBluesteinTest.java new file mode 100644 index 00000000..9127e16c --- /dev/null +++ b/src/test/java/com/maths/FFTBluesteinTest.java @@ -0,0 +1,79 @@ +package com.maths; + +import org.junit.jupiter.api.Test; +import java.math.BigDecimal; +import java.math.RoundingMode; +import java.util.ArrayList; +import static org.junit.jupiter.api.Assertions.*; + +/** + * Class for testing the Bluestein's FFT algorithm. + * You can also use Matlab for bigger signals and compare the results by modifying the code below. + * + * @author Ioannis Karavitsis + * @version 1.0 + * */ +class FFTBluesteinTest +{ + /** + * Function to round a number with double precision to n decimal places. + * + * More info: + * https://www.baeldung.com/java-round-decimal-number + * + * @param num The number to be rounded. + * @param places The number of decimal places up to which the number will be rounded. + * @return The rounded number. + * */ + private double round(double num, int places) + { + if (places < 0) throw new IllegalArgumentException(); + + BigDecimal bd = new BigDecimal(Double.toString(num)); + bd = bd.setScale(places, RoundingMode.HALF_UP); + return bd.doubleValue(); + } + + /** + * Test function for fftBluestein() function. + * */ + @Test + public void testFFT() + { + ArrayList x = new ArrayList<>(); + for (int i = 1; i <= 10; i++) //x = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 } + x.add(new FFT.Complex(i, 0)); + + /* Matlab results */ + double[] real = { 55, -5, -5, -5, -5, -5, -5, -5, -5, -5 }; + double[] imaginary = { 0, 15.388417685876266, 6.881909602355867, 3.632712640026803, 1.624598481164532, 0, -1.624598481164532, -3.632712640026803, -6.881909602355867, -15.388417685876266 }; + + FFTBluestein.fftBluestein(x,false); + + /* Print the results (Optional) */ + /*for(int i = 0; i < x.size(); i++) + System.out.println(x.get(i).getReal() + " " + x.get(i).getImaginary());*/ + + /* Compare the results (real and imaginary part) with those from Matlab (command: fft(x)). + I rounded them to 12 decimal places, due to the different calculations and precision that Matlab and this Java program use. */ + for(int i = 0; i < x.size(); i++) + { + assertEquals(round(real[i], 0), round(x.get(i).getReal(), 0)); + assertEquals(round(imaginary[i], 12), round(x.get(i).getImaginary(), 12)); + } + + /* Find the IFFT of the FFT of signal x */ + FFTBluestein.fftBluestein(x, true); + + /* Print the results (Optional) */ + /*for(int i = 0; i < x.size(); i++) + System.out.println(x.get(i).getReal() + " " + x.get(i).getImaginary());*/ + + /* Check if it's equal to signal x */ + for(int i = 0; i < 10; i++) + { + assertEquals(i + 1, round(x.get(i).getReal(), 0)); + assertEquals(0, round(x.get(i).getImaginary(), 0)); //All imaginary parts must be 0. + } + } +} \ No newline at end of file diff --git a/src/test/java/com/maths/FFTTest.java b/src/test/java/com/maths/FFTTest.java new file mode 100644 index 00000000..7b3fffe7 --- /dev/null +++ b/src/test/java/com/maths/FFTTest.java @@ -0,0 +1,79 @@ +package com.maths; + +import org.junit.jupiter.api.Test; +import java.math.BigDecimal; +import java.math.RoundingMode; +import java.util.ArrayList; +import static org.junit.jupiter.api.Assertions.*; + +/** + * Class for testing the FFT Cooley-Tukey algorithm. + * You can also use Matlab for bigger signals and compare the results by modifying the code below. + * + * @author Ioannis Karavitsis + * @version 1.0 + * */ +class FFTTest +{ + /** + * Function to round a number with double precision to n decimal places. + * + * More info: + * https://www.baeldung.com/java-round-decimal-number + * + * @param num The number to be rounded. + * @param places The number of decimal places up to which the number will be rounded. + * @return The rounded number. + * */ + private double round(double num, int places) + { + if (places < 0) throw new IllegalArgumentException(); + + BigDecimal bd = new BigDecimal(Double.toString(num)); + bd = bd.setScale(places, RoundingMode.HALF_UP); + return bd.doubleValue(); + } + + /** + * Test function for fft() function. + * */ + @Test + public void testFFT() + { + ArrayList x = new ArrayList<>(); + for (int i = 1; i <= 10; i++) //x = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 } + x.add(new FFT.Complex(i, 0)); + + /* Matlab results */ + double[] real = { 55, -26.375866509656959, 12.071067811865476, -9.446748728072674, 5, -0.896397022434947, -2.071067811865476, 4.719012260164577, -5, 4.719012260164577, -2.071067811865476, -0.896397022434947, 5, -9.446748728072674, 12.071067811865476, -26.37586650965695 }; + double[] imaginary = { 0, -21.309863136978343, 2.585786437626905, 1.755766511785423, -6, 5.897902135516374, -5.414213562373095, 2.832272486752609, 0, -2.832272486752609, 5.414213562373095, -5.897902135516374, 6, -1.755766511785423, -2.585786437626905, 21.309863136978343 }; + + FFT.fft(x,false); + + /* Print the results (Optional) */ + /*for(int i = 0; i < x.size(); i++) + System.out.println(x.get(i).getReal() + " " + x.get(i).getImaginary());*/ + + /* Compare the results (real and imaginary part) with those from Matlab (command: fft(x, 16)). + I rounded them to 13 decimal places, due to the different calculations and precision that Matlab and this Java program use. */ + for(int i = 0; i < x.size(); i++) + { + assertEquals(round(real[i], 13), round(x.get(i).getReal(), 13)); + assertEquals(round(imaginary[i], 13), round(x.get(i).getImaginary(), 13)); + } + + /* Find the IFFT of the FFT of signal x */ + FFT.fft(x,true); + + /* Print the results (Optional) */ + /*for(int i = 0; i < x.size(); i++) + System.out.println(x.get(i).getReal() + " " + x.get(i).getImaginary());*/ + + /* Check if it's equal to signal x */ + for(int i = 0; i < 10; i++) //The rest values are 0, so we don't have to check them. + { + assertEquals(i + 1, round(x.get(i).getReal(), 0)); + assertEquals(0, round(x.get(i).getImaginary(), 0)); //All imaginary parts must be 0. + } + } +} \ No newline at end of file