Formatted with Google Java Formatter
This commit is contained in:
parent
4264e0f045
commit
ca2e20707f
@ -7,40 +7,38 @@ import java.util.ArrayList;
|
|||||||
*
|
*
|
||||||
* @author Ioannis Karavitsis
|
* @author Ioannis Karavitsis
|
||||||
* @version 1.0
|
* @version 1.0
|
||||||
* */
|
*/
|
||||||
public class CircularConvolutionFFT
|
public class CircularConvolutionFFT {
|
||||||
{
|
|
||||||
/**
|
/**
|
||||||
* This method pads the signal with zeros until it reaches the new size.
|
* This method pads the signal with zeros until it reaches the new size.
|
||||||
*
|
*
|
||||||
* @param x The signal to be padded.
|
* @param x The signal to be padded.
|
||||||
* @param newSize The new size of the signal.
|
* @param newSize The new size of the signal.
|
||||||
* */
|
*/
|
||||||
private static void padding(ArrayList<FFT.Complex> x, int newSize)
|
private static void padding(ArrayList<FFT.Complex> x, int newSize) {
|
||||||
{
|
if (x.size() < newSize) {
|
||||||
if(x.size() < newSize)
|
|
||||||
{
|
|
||||||
int diff = newSize - x.size();
|
int diff = newSize - x.size();
|
||||||
for(int i = 0; i < diff; i++)
|
for (int i = 0; i < diff; i++) x.add(new FFT.Complex());
|
||||||
x.add(new FFT.Complex());
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Discrete circular convolution function. It uses the convolution theorem for discrete signals: convolved = IDFT(DFT(a)*DFT(b)).
|
* Discrete circular convolution function. It uses the convolution theorem for discrete signals:
|
||||||
* Then we use the FFT algorithm for faster calculations of the two DFTs and the final IDFT.
|
* 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:
|
* <p>More info: https://en.wikipedia.org/wiki/Convolution_theorem
|
||||||
* https://en.wikipedia.org/wiki/Convolution_theorem
|
|
||||||
*
|
*
|
||||||
* @param a The first signal.
|
* @param a The first signal.
|
||||||
* @param b The other signal.
|
* @param b The other signal.
|
||||||
* @return The convolved signal.
|
* @return The convolved signal.
|
||||||
* */
|
*/
|
||||||
public static ArrayList<FFT.Complex> fftCircularConvolution(ArrayList<FFT.Complex> a, ArrayList<FFT.Complex> b)
|
public static ArrayList<FFT.Complex> fftCircularConvolution(
|
||||||
{
|
ArrayList<FFT.Complex> a, ArrayList<FFT.Complex> b) {
|
||||||
int convolvedSize = Math.max(a.size(), b.size()); //The two signals must have the same size equal to the bigger one
|
int convolvedSize =
|
||||||
padding(a, convolvedSize); //Zero padding the smaller signal
|
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);
|
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 */
|
/* 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 */
|
||||||
@ -48,10 +46,9 @@ public class CircularConvolutionFFT
|
|||||||
FFTBluestein.fftBluestein(b, false);
|
FFTBluestein.fftBluestein(b, false);
|
||||||
ArrayList<FFT.Complex> convolved = new ArrayList<>();
|
ArrayList<FFT.Complex> convolved = new ArrayList<>();
|
||||||
|
|
||||||
for(int i = 0; i < a.size(); i++)
|
for (int i = 0; i < a.size(); i++) convolved.add(a.get(i).multiply(b.get(i))); // FFT(a)*FFT(b)
|
||||||
convolved.add(a.get(i).multiply(b.get(i))); //FFT(a)*FFT(b)
|
|
||||||
|
|
||||||
FFTBluestein.fftBluestein(convolved, true); //IFFT
|
FFTBluestein.fftBluestein(convolved, true); // IFFT
|
||||||
|
|
||||||
return convolved;
|
return convolved;
|
||||||
}
|
}
|
||||||
|
@ -5,19 +5,17 @@ package com.maths;
|
|||||||
*
|
*
|
||||||
* @author Ioannis Karavitsis
|
* @author Ioannis Karavitsis
|
||||||
* @version 1.0
|
* @version 1.0
|
||||||
* */
|
*/
|
||||||
public class Convolution
|
public class Convolution {
|
||||||
{
|
|
||||||
/**
|
/**
|
||||||
* Discrete linear convolution function. Both input signals and the output signal must start from 0.
|
* Discrete linear convolution function. Both input signals and the output signal must start from
|
||||||
* If you have a signal that has values before 0 then shift it to start from 0.
|
* 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 A The first discrete signal
|
||||||
* @param B The second discrete signal
|
* @param B The second discrete signal
|
||||||
* @return The convolved signal
|
* @return The convolved signal
|
||||||
* */
|
*/
|
||||||
public static double[] convolution(double[] A, double[] B)
|
public static double[] convolution(double[] A, double[] B) {
|
||||||
{
|
|
||||||
double[] convolved = new double[A.length + B.length - 1];
|
double[] convolved = new double[A.length + B.length - 1];
|
||||||
|
|
||||||
/*
|
/*
|
||||||
@ -30,13 +28,11 @@ public class Convolution
|
|||||||
It's obvious that: 0 <= k <= A.length , 0 <= i <= A.length + B.length - 2 and 0 <= i-k <= B.length - 1
|
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.
|
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++)
|
for (int i = 0; i < convolved.length; i++) {
|
||||||
{
|
|
||||||
convolved[i] = 0;
|
convolved[i] = 0;
|
||||||
int k = Math.max(i - B.length + 1, 0);
|
int k = Math.max(i - B.length + 1, 0);
|
||||||
|
|
||||||
while(k < i + 1 && k < A.length)
|
while (k < i + 1 && k < A.length) {
|
||||||
{
|
|
||||||
convolved[i] += A[k] * B[i - k];
|
convolved[i] += A[k] * B[i - k];
|
||||||
k++;
|
k++;
|
||||||
}
|
}
|
||||||
|
@ -7,42 +7,39 @@ import java.util.ArrayList;
|
|||||||
*
|
*
|
||||||
* @author Ioannis Karavitsis
|
* @author Ioannis Karavitsis
|
||||||
* @version 1.0
|
* @version 1.0
|
||||||
* */
|
*/
|
||||||
public class ConvolutionFFT
|
public class ConvolutionFFT {
|
||||||
{
|
|
||||||
/**
|
/**
|
||||||
* This method pads the signal with zeros until it reaches the new size.
|
* This method pads the signal with zeros until it reaches the new size.
|
||||||
*
|
*
|
||||||
* @param x The signal to be padded.
|
* @param x The signal to be padded.
|
||||||
* @param newSize The new size of the signal.
|
* @param newSize The new size of the signal.
|
||||||
* */
|
*/
|
||||||
private static void padding(ArrayList<FFT.Complex> x, int newSize)
|
private static void padding(ArrayList<FFT.Complex> x, int newSize) {
|
||||||
{
|
if (x.size() < newSize) {
|
||||||
if(x.size() < newSize)
|
|
||||||
{
|
|
||||||
int diff = newSize - x.size();
|
int diff = newSize - x.size();
|
||||||
for(int i = 0; i < diff; i++)
|
for (int i = 0; i < diff; i++) x.add(new FFT.Complex());
|
||||||
x.add(new FFT.Complex());
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Discrete linear convolution function. It uses the convolution theorem for discrete signals convolved: = IDFT(DFT(a)*DFT(b)).
|
* Discrete linear convolution function. It uses the convolution theorem for discrete signals
|
||||||
* 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).
|
* convolved: = IDFT(DFT(a)*DFT(b)). This is true for circular convolution. In order to get the
|
||||||
* Then we use the FFT algorithm for faster calculations of the two DFTs and the final IDFT.
|
* 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:
|
* <p>More info: https://en.wikipedia.org/wiki/Convolution_theorem
|
||||||
* https://en.wikipedia.org/wiki/Convolution_theorem
|
|
||||||
* https://ccrma.stanford.edu/~jos/ReviewFourier/FFT_Convolution.html
|
* https://ccrma.stanford.edu/~jos/ReviewFourier/FFT_Convolution.html
|
||||||
*
|
*
|
||||||
* @param a The first signal.
|
* @param a The first signal.
|
||||||
* @param b The other signal.
|
* @param b The other signal.
|
||||||
* @return The convolved signal.
|
* @return The convolved signal.
|
||||||
* */
|
*/
|
||||||
public static ArrayList<FFT.Complex> convolutionFFT(ArrayList<FFT.Complex> a, ArrayList<FFT.Complex> b)
|
public static ArrayList<FFT.Complex> convolutionFFT(
|
||||||
{
|
ArrayList<FFT.Complex> a, ArrayList<FFT.Complex> b) {
|
||||||
int convolvedSize = a.size() + b.size() - 1; //The size of the convolved signal
|
int convolvedSize = a.size() + b.size() - 1; // The size of the convolved signal
|
||||||
padding(a, convolvedSize); //Zero padding both signals
|
padding(a, convolvedSize); // Zero padding both signals
|
||||||
padding(b, convolvedSize);
|
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) */
|
/* 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) */
|
||||||
@ -50,11 +47,13 @@ public class ConvolutionFFT
|
|||||||
FFT.fft(b, false);
|
FFT.fft(b, false);
|
||||||
ArrayList<FFT.Complex> convolved = new ArrayList<>();
|
ArrayList<FFT.Complex> convolved = new ArrayList<>();
|
||||||
|
|
||||||
for(int i = 0; i < a.size(); i++)
|
for (int i = 0; i < a.size(); i++) convolved.add(a.get(i).multiply(b.get(i))); // FFT(a)*FFT(b)
|
||||||
convolved.add(a.get(i).multiply(b.get(i))); //FFT(a)*FFT(b)
|
|
||||||
|
|
||||||
FFT.fft(convolved, true); //IFFT
|
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.
|
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;
|
return convolved;
|
||||||
}
|
}
|
||||||
|
164
Maths/FFT.java
164
Maths/FFT.java
@ -4,29 +4,23 @@ import java.util.ArrayList;
|
|||||||
import java.util.Collections;
|
import java.util.Collections;
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Class for calculating the Fast Fourier Transform (FFT) of a discrete signal using the Cooley-Tukey algorithm.
|
* Class for calculating the Fast Fourier Transform (FFT) of a discrete signal using the
|
||||||
|
* Cooley-Tukey algorithm.
|
||||||
*
|
*
|
||||||
* @author Ioannis Karavitsis
|
* @author Ioannis Karavitsis
|
||||||
* @version 1.0
|
* @version 1.0
|
||||||
* */
|
*/
|
||||||
public class FFT
|
public class FFT {
|
||||||
{
|
|
||||||
/**
|
/**
|
||||||
* This class represents a complex number and has methods for basic operations.
|
* This class represents a complex number and has methods for basic operations.
|
||||||
*
|
*
|
||||||
* More info:
|
* <p>More info: https://introcs.cs.princeton.edu/java/32class/Complex.java.html
|
||||||
* https://introcs.cs.princeton.edu/java/32class/Complex.java.html
|
*/
|
||||||
* */
|
static class Complex {
|
||||||
static class Complex
|
|
||||||
{
|
|
||||||
private double real, img;
|
private double real, img;
|
||||||
|
|
||||||
/**
|
/** Default Constructor. Creates the complex number 0. */
|
||||||
* Default Constructor.
|
public Complex() {
|
||||||
* Creates the complex number 0.
|
|
||||||
* */
|
|
||||||
public Complex()
|
|
||||||
{
|
|
||||||
real = 0;
|
real = 0;
|
||||||
img = 0;
|
img = 0;
|
||||||
}
|
}
|
||||||
@ -36,9 +30,8 @@ public class FFT
|
|||||||
*
|
*
|
||||||
* @param r The real part of the number.
|
* @param r The real part of the number.
|
||||||
* @param i The imaginary part of the number.
|
* @param i The imaginary part of the number.
|
||||||
* */
|
*/
|
||||||
public Complex(double r, double i)
|
public Complex(double r, double i) {
|
||||||
{
|
|
||||||
real = r;
|
real = r;
|
||||||
img = i;
|
img = i;
|
||||||
}
|
}
|
||||||
@ -47,9 +40,8 @@ public class FFT
|
|||||||
* Returns the real part of the complex number.
|
* Returns the real part of the complex number.
|
||||||
*
|
*
|
||||||
* @return The real part of the complex number.
|
* @return The real part of the complex number.
|
||||||
* */
|
*/
|
||||||
public double getReal()
|
public double getReal() {
|
||||||
{
|
|
||||||
return real;
|
return real;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -57,9 +49,8 @@ public class FFT
|
|||||||
* Returns the imaginary part of the complex number.
|
* Returns the imaginary part of the complex number.
|
||||||
*
|
*
|
||||||
* @return The imaginary part of the complex number.
|
* @return The imaginary part of the complex number.
|
||||||
* */
|
*/
|
||||||
public double getImaginary()
|
public double getImaginary() {
|
||||||
{
|
|
||||||
return img;
|
return img;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -68,9 +59,8 @@ public class FFT
|
|||||||
*
|
*
|
||||||
* @param z The number to be added.
|
* @param z The number to be added.
|
||||||
* @return The sum.
|
* @return The sum.
|
||||||
* */
|
*/
|
||||||
public Complex add(Complex z)
|
public Complex add(Complex z) {
|
||||||
{
|
|
||||||
Complex temp = new Complex();
|
Complex temp = new Complex();
|
||||||
temp.real = this.real + z.real;
|
temp.real = this.real + z.real;
|
||||||
temp.img = this.img + z.img;
|
temp.img = this.img + z.img;
|
||||||
@ -82,9 +72,8 @@ public class FFT
|
|||||||
*
|
*
|
||||||
* @param z The number to be subtracted.
|
* @param z The number to be subtracted.
|
||||||
* @return The difference.
|
* @return The difference.
|
||||||
* */
|
*/
|
||||||
public Complex subtract(Complex z)
|
public Complex subtract(Complex z) {
|
||||||
{
|
|
||||||
Complex temp = new Complex();
|
Complex temp = new Complex();
|
||||||
temp.real = this.real - z.real;
|
temp.real = this.real - z.real;
|
||||||
temp.img = this.img - z.img;
|
temp.img = this.img - z.img;
|
||||||
@ -96,12 +85,11 @@ public class FFT
|
|||||||
*
|
*
|
||||||
* @param z The number to be multiplied.
|
* @param z The number to be multiplied.
|
||||||
* @return The product.
|
* @return The product.
|
||||||
* */
|
*/
|
||||||
public Complex multiply(Complex z)
|
public Complex multiply(Complex z) {
|
||||||
{
|
|
||||||
Complex temp = new Complex();
|
Complex temp = new Complex();
|
||||||
temp.real = this.real*z.real - this.img*z.img;
|
temp.real = this.real * z.real - this.img * z.img;
|
||||||
temp.img = this.real*z.img + this.img*z.real;
|
temp.img = this.real * z.img + this.img * z.real;
|
||||||
return temp;
|
return temp;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -110,9 +98,8 @@ public class FFT
|
|||||||
*
|
*
|
||||||
* @param n The real number to be multiplied.
|
* @param n The real number to be multiplied.
|
||||||
* @return The product.
|
* @return The product.
|
||||||
* */
|
*/
|
||||||
public Complex multiply(double n)
|
public Complex multiply(double n) {
|
||||||
{
|
|
||||||
Complex temp = new Complex();
|
Complex temp = new Complex();
|
||||||
temp.real = this.real * n;
|
temp.real = this.real * n;
|
||||||
temp.img = this.img * n;
|
temp.img = this.img * n;
|
||||||
@ -123,9 +110,8 @@ public class FFT
|
|||||||
* Finds the conjugate of this complex number.
|
* Finds the conjugate of this complex number.
|
||||||
*
|
*
|
||||||
* @return The conjugate.
|
* @return The conjugate.
|
||||||
* */
|
*/
|
||||||
public Complex conjugate()
|
public Complex conjugate() {
|
||||||
{
|
|
||||||
Complex temp = new Complex();
|
Complex temp = new Complex();
|
||||||
temp.real = this.real;
|
temp.real = this.real;
|
||||||
temp.img = -this.img;
|
temp.img = -this.img;
|
||||||
@ -136,9 +122,8 @@ public class FFT
|
|||||||
* Finds the magnitude of the complex number.
|
* Finds the magnitude of the complex number.
|
||||||
*
|
*
|
||||||
* @return The magnitude.
|
* @return The magnitude.
|
||||||
* */
|
*/
|
||||||
public double abs()
|
public double abs() {
|
||||||
{
|
|
||||||
return Math.hypot(this.real, this.img);
|
return Math.hypot(this.real, this.img);
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -147,12 +132,11 @@ public class FFT
|
|||||||
*
|
*
|
||||||
* @param z The divisor.
|
* @param z The divisor.
|
||||||
* @return The quotient.
|
* @return The quotient.
|
||||||
* */
|
*/
|
||||||
public Complex divide(Complex z)
|
public Complex divide(Complex z) {
|
||||||
{
|
|
||||||
Complex temp = new Complex();
|
Complex temp = new Complex();
|
||||||
temp.real = (this.real*z.real + this.img*z.img) / (z.abs()*z.abs());
|
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());
|
temp.img = (this.img * z.real - this.real * z.img) / (z.abs() * z.abs());
|
||||||
return temp;
|
return temp;
|
||||||
}
|
}
|
||||||
|
|
||||||
@ -161,9 +145,8 @@ public class FFT
|
|||||||
*
|
*
|
||||||
* @param n The divisor which is a real number.
|
* @param n The divisor which is a real number.
|
||||||
* @return The quotient.
|
* @return The quotient.
|
||||||
* */
|
*/
|
||||||
public Complex divide(double n)
|
public Complex divide(double n) {
|
||||||
{
|
|
||||||
Complex temp = new Complex();
|
Complex temp = new Complex();
|
||||||
temp.real = this.real / n;
|
temp.real = this.real / n;
|
||||||
temp.img = this.img / n;
|
temp.img = this.img / n;
|
||||||
@ -172,64 +155,55 @@ public class FFT
|
|||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Iterative In-Place Radix-2 Cooley-Tukey Fast Fourier Transform Algorithm with Bit-Reversal.
|
* Iterative In-Place Radix-2 Cooley-Tukey Fast Fourier Transform Algorithm with Bit-Reversal. The
|
||||||
* 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.
|
* 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:
|
* <p>More info: https://www.algorithm-archive.org/contents/cooley_tukey/cooley_tukey.html
|
||||||
* https://www.algorithm-archive.org/contents/cooley_tukey/cooley_tukey.html
|
|
||||||
* https://www.geeksforgeeks.org/iterative-fast-fourier-transformation-polynomial-multiplication/
|
* https://www.geeksforgeeks.org/iterative-fast-fourier-transformation-polynomial-multiplication/
|
||||||
* https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm
|
* https://en.wikipedia.org/wiki/Cooley%E2%80%93Tukey_FFT_algorithm
|
||||||
* https://cp-algorithms.com/algebra/fft.html
|
* 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 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.
|
* @param inverse True if you want to find the inverse FFT.
|
||||||
* */
|
*/
|
||||||
public static void fft(ArrayList<Complex> x, boolean inverse)
|
public static void fft(ArrayList<Complex> x, boolean inverse) {
|
||||||
{
|
|
||||||
/* Pad the signal with zeros if necessary */
|
/* Pad the signal with zeros if necessary */
|
||||||
paddingPowerOfTwo(x);
|
paddingPowerOfTwo(x);
|
||||||
int N = x.size();
|
int N = x.size();
|
||||||
|
|
||||||
/* Find the log2(N) */
|
/* Find the log2(N) */
|
||||||
int log2N = 0;
|
int log2N = 0;
|
||||||
while((1 << log2N) < N)
|
while ((1 << log2N) < N) log2N++;
|
||||||
log2N++;
|
|
||||||
|
|
||||||
/* Swap the values of the signal with bit-reversal method */
|
/* Swap the values of the signal with bit-reversal method */
|
||||||
int reverse;
|
int reverse;
|
||||||
for(int i = 0; i < N; i++)
|
for (int i = 0; i < N; i++) {
|
||||||
{
|
|
||||||
reverse = reverseBits(i, log2N);
|
reverse = reverseBits(i, log2N);
|
||||||
if(i < reverse)
|
if (i < reverse) Collections.swap(x, i, reverse);
|
||||||
Collections.swap(x, i, reverse);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
int direction = inverse ? -1 : 1;
|
int direction = inverse ? -1 : 1;
|
||||||
|
|
||||||
/* Main loop of the algorithm */
|
/* Main loop of the algorithm */
|
||||||
for(int len = 2; len <= N; len *= 2)
|
for (int len = 2; len <= N; len *= 2) {
|
||||||
{
|
|
||||||
double angle = -2 * Math.PI / len * direction;
|
double angle = -2 * Math.PI / len * direction;
|
||||||
Complex wlen = new Complex(Math.cos(angle), Math.sin(angle));
|
Complex wlen = new Complex(Math.cos(angle), Math.sin(angle));
|
||||||
for(int i = 0; i < N; i += len)
|
for (int i = 0; i < N; i += len) {
|
||||||
{
|
|
||||||
Complex w = new Complex(1, 0);
|
Complex w = new Complex(1, 0);
|
||||||
for(int j = 0; j < len / 2; j++)
|
for (int j = 0; j < len / 2; j++) {
|
||||||
{
|
|
||||||
Complex u = x.get(i + j);
|
Complex u = x.get(i + j);
|
||||||
Complex v = w.multiply(x.get(i + j + len/2));
|
Complex v = w.multiply(x.get(i + j + len / 2));
|
||||||
x.set(i + j, u.add(v));
|
x.set(i + j, u.add(v));
|
||||||
x.set(i + j + len/2, u.subtract(v));
|
x.set(i + j + len / 2, u.subtract(v));
|
||||||
w = w.multiply(wlen);
|
w = w.multiply(wlen);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
/* Divide by N if we want the inverse FFT */
|
/* Divide by N if we want the inverse FFT */
|
||||||
if(inverse)
|
if (inverse) {
|
||||||
{
|
for (int i = 0; i < x.size(); i++) {
|
||||||
for (int i = 0; i < x.size(); i++)
|
|
||||||
{
|
|
||||||
Complex z = x.get(i);
|
Complex z = x.get(i);
|
||||||
x.set(i, z.divide(N));
|
x.set(i, z.divide(N));
|
||||||
}
|
}
|
||||||
@ -237,45 +211,35 @@ public class FFT
|
|||||||
}
|
}
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* This function reverses the bits of a number.
|
* This function reverses the bits of a number. It is used in Cooley-Tukey FFT algorithm.
|
||||||
* It is used in Cooley-Tukey FFT algorithm.
|
|
||||||
*
|
*
|
||||||
* E.g.
|
* <p>E.g. num = 13 = 00001101 in binary log2N = 8 Then reversed = 176 = 10110000 in binary
|
||||||
* num = 13 = 00001101 in binary
|
|
||||||
* log2N = 8
|
|
||||||
* Then reversed = 176 = 10110000 in binary
|
|
||||||
*
|
*
|
||||||
* More info:
|
* <p>More info: https://cp-algorithms.com/algebra/fft.html
|
||||||
* https://cp-algorithms.com/algebra/fft.html
|
|
||||||
* https://www.geeksforgeeks.org/write-an-efficient-c-program-to-reverse-bits-of-a-number/
|
* 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 num The integer you want to reverse its bits.
|
||||||
* @param log2N The number of bits you want to reverse.
|
* @param log2N The number of bits you want to reverse.
|
||||||
* @return The reversed number
|
* @return The reversed number
|
||||||
* */
|
*/
|
||||||
private static int reverseBits(int num, int log2N)
|
private static int reverseBits(int num, int log2N) {
|
||||||
{
|
|
||||||
int reversed = 0;
|
int reversed = 0;
|
||||||
for(int i = 0; i < log2N; i++)
|
for (int i = 0; i < log2N; i++) {
|
||||||
{
|
if ((num & (1 << i)) != 0) reversed |= 1 << (log2N - 1 - i);
|
||||||
if((num & (1 << i)) != 0)
|
|
||||||
reversed |= 1 << (log2N - 1 - i);
|
|
||||||
}
|
}
|
||||||
return reversed;
|
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.
|
* 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.
|
* @param x The ArrayList to be padded.
|
||||||
* */
|
*/
|
||||||
private static void paddingPowerOfTwo(ArrayList<Complex> x)
|
private static void paddingPowerOfTwo(ArrayList<Complex> x) {
|
||||||
{
|
|
||||||
int n = 1;
|
int n = 1;
|
||||||
int oldSize = x.size();
|
int oldSize = x.size();
|
||||||
while(n < oldSize)
|
while (n < oldSize) n *= 2;
|
||||||
n *= 2;
|
for (int i = 0; i < n - oldSize; i++) x.add(new Complex());
|
||||||
for(int i = 0; i < n - oldSize; i++)
|
|
||||||
x.add(new Complex());
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -3,45 +3,40 @@ package com.maths;
|
|||||||
import java.util.ArrayList;
|
import java.util.ArrayList;
|
||||||
|
|
||||||
/**
|
/**
|
||||||
* Class for calculating the Fast Fourier Transform (FFT) of a discrete signal using the Bluestein's algorithm.
|
* Class for calculating the Fast Fourier Transform (FFT) of a discrete signal using the Bluestein's
|
||||||
|
* algorithm.
|
||||||
*
|
*
|
||||||
* @author Ioannis Karavitsis
|
* @author Ioannis Karavitsis
|
||||||
* @version 1.0
|
* @version 1.0
|
||||||
* */
|
*/
|
||||||
public class FFTBluestein
|
public class FFTBluestein {
|
||||||
{
|
|
||||||
/**
|
/**
|
||||||
* Bluestein's FFT Algorithm.
|
* Bluestein's FFT Algorithm.
|
||||||
*
|
*
|
||||||
* More info:
|
* <p>More info: https://en.wikipedia.org/wiki/Chirp_Z-transform#Bluestein.27s_algorithm
|
||||||
* 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
|
* 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 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.
|
* @param inverse True if you want to find the inverse FFT.
|
||||||
* */
|
*/
|
||||||
public static void fftBluestein(ArrayList<FFT.Complex> x, boolean inverse)
|
public static void fftBluestein(ArrayList<FFT.Complex> x, boolean inverse) {
|
||||||
{
|
|
||||||
int N = x.size();
|
int N = x.size();
|
||||||
int bnSize = 2*N - 1;
|
int bnSize = 2 * N - 1;
|
||||||
int direction = inverse ? -1 : 1;
|
int direction = inverse ? -1 : 1;
|
||||||
ArrayList<FFT.Complex> an = new ArrayList<>();
|
ArrayList<FFT.Complex> an = new ArrayList<>();
|
||||||
ArrayList<FFT.Complex> bn = new ArrayList<>();
|
ArrayList<FFT.Complex> bn = new ArrayList<>();
|
||||||
|
|
||||||
/* Initialization of the b(n) sequence (see Wikipedia's article above for the symbols used)*/
|
/* Initialization of the b(n) sequence (see Wikipedia's article above for the symbols used)*/
|
||||||
for(int i = 0; i < bnSize; i++)
|
for (int i = 0; i < bnSize; i++) bn.add(new FFT.Complex());
|
||||||
bn.add(new FFT.Complex());
|
|
||||||
|
|
||||||
for(int i = 0; i < N; i++)
|
for (int i = 0; i < N; i++) {
|
||||||
{
|
|
||||||
double angle = (i - N + 1) * (i - N + 1) * Math.PI / N * direction;
|
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(i, new FFT.Complex(Math.cos(angle), Math.sin(angle)));
|
||||||
bn.set(bnSize - i - 1, 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 */
|
/* Initialization of the a(n) sequence */
|
||||||
for(int i = 0; i < N; i++)
|
for (int i = 0; i < N; i++) {
|
||||||
{
|
|
||||||
double angle = -i * i * Math.PI / N * direction;
|
double angle = -i * i * Math.PI / N * direction;
|
||||||
an.add(x.get(i).multiply(new FFT.Complex(Math.cos(angle), Math.sin(angle))));
|
an.add(x.get(i).multiply(new FFT.Complex(Math.cos(angle), Math.sin(angle))));
|
||||||
}
|
}
|
||||||
@ -49,18 +44,15 @@ public class FFTBluestein
|
|||||||
ArrayList<FFT.Complex> convolution = ConvolutionFFT.convolutionFFT(an, bn);
|
ArrayList<FFT.Complex> convolution = ConvolutionFFT.convolutionFFT(an, bn);
|
||||||
|
|
||||||
/* The final multiplication of the convolution with the b*(k) factor */
|
/* The final multiplication of the convolution with the b*(k) factor */
|
||||||
for(int i = 0; i < N; i++)
|
for (int i = 0; i < N; i++) {
|
||||||
{
|
|
||||||
double angle = -1 * i * i * Math.PI / N * direction;
|
double angle = -1 * i * i * Math.PI / N * direction;
|
||||||
FFT.Complex bk = new FFT.Complex(Math.cos(angle), Math.sin(angle));
|
FFT.Complex bk = new FFT.Complex(Math.cos(angle), Math.sin(angle));
|
||||||
x.set(i, bk.multiply(convolution.get(i + N - 1)));
|
x.set(i, bk.multiply(convolution.get(i + N - 1)));
|
||||||
}
|
}
|
||||||
|
|
||||||
/* Divide by N if we want the inverse FFT */
|
/* Divide by N if we want the inverse FFT */
|
||||||
if(inverse)
|
if (inverse) {
|
||||||
{
|
for (int i = 0; i < N; i++) {
|
||||||
for (int i = 0; i < N; i++)
|
|
||||||
{
|
|
||||||
FFT.Complex z = x.get(i);
|
FFT.Complex z = x.get(i);
|
||||||
x.set(i, z.divide(N));
|
x.set(i, z.divide(N));
|
||||||
}
|
}
|
||||||
|
Loading…
Reference in New Issue
Block a user