Fast Fourier transform(고속 퓨리에 변환)
고속 퓨리에 변환FFT는 이산 퓨리에 변환DFT와 그 역변환을 빠르게 수행하는 효율적인 알고리즘 입니다.
이 복소수complex 일 때, DFT는 다음과 같습니다.
위 식에 따르면 O(n^2)의 연산이 필요하지만, FFT를 이용하면 O(n log n)의 연산만으로 가능합니다.
가장 일반적으로 사용되는 FFT 알고리즘은 쿨리-튜키 알고리즘입니다.
참고로, DFT(Discrete Fourier Transform)는 이산화된 시간 영역의 데이터를 이산화된 주파수 영역으로 변환해주는 알고리즘입니다. 일반적인 디지털 신호를 디지털 주파수로 바꾸어줍니다.
왼쪽으로 보면 시간축을 기준으로 하는 주기함수들의 합이 붉은 선의 모양으로 나타나게 됩니다.
하지만 오른쪽으로 보면 주파수를 기준으로 분리된 각각의 주파수 영역과 크기를 볼 수 있습니다.
디지털 신호를 안드로이드에서 FFT 하기 위하여, 3가지의 FFT 라이브러리와 코드를 사용해보고 비교해보았습니다.
Columbia University FFT class
제일 먼저 찾은 것은 Columbia University 출처의 FFT 코드 입니다.
public class FFT {
int n, m;
// Lookup tables. Only need to recompute when size of FFT changes.
double[] cos;
double[] sin;
double[] window;
public FFT(int n) {
this.n = n;
this.m = (int)(Math.log(n) / Math.log(2));
// Make sure n is a power of 2
if(n != (1<<m))
throw new RuntimeException("FFT length must be power of 2");
// precompute tables
cos = new double[n/2];
sin = new double[n/2];
for(int i=0; i<n/2; i++) {
cos[i] = Math.cos(-2*Math.PI*i/n);
sin[i] = Math.sin(-2*Math.PI*i/n);
}
}
public void fft(double[] x, double[] y)
{
int i,j,k,n1,n2,a;
double c,s,e,t1,t2;
// Bit-reverse
j = 0;
n2 = n/2;
for (i=1; i < n - 1; i++) {
n1 = n2;
while ( j >= n1 ) {
j = j - n1;
n1 = n1/2;
}
j = j + n1;
if (i < j) {
t1 = x[i];
x[i] = x[j];
x[j] = t1;
t1 = y[i];
y[i] = y[j];
y[j] = t1;
}
}
// FFT
n1 = 0;
n2 = 1;
for (i=0; i < m; i++) {
n1 = n2;
n2 = n2 + n2;
a = 0;
for (j=0; j < n1; j++) {
c = cos[a];
s = sin[a];
a += 1 << (m-i-1);
for (k=j; k < n; k=k+n2) {
t1 = c*x[k+n1] - s*y[k+n1];
t2 = s*x[k+n1] + c*y[k+n1];
x[k+n1] = x[k] - t1;
y[k+n1] = y[k] - t2;
x[k] = x[k] + t1;
y[k] = y[k] + t2;
}
}
}
}
}
코드 활용 방법
위의 class 또는 함수를 사용하면, input을 FFT output으로 바꿀 수 있습니다.
Input
- n : the number of data points ( input array의 사이즈, 2의 제곱이여야 함.)
- x : transform될 데이터의 실수 부분
- y : transform될 데이터의 허수 부분
Output
- x : FFT output의 실수 부분
- y : FFT output의 허수 부분
일반적인 FFT graph를 그리기 위해서는 magnitude를 계산하여야 합니다.
double[] mag = new double[N/2];
for(int i=0;i<N/2;i++) {
mag[i] = Math.sqrt(Math.pow(x[i],2)+Math.pow(y[i],2));
}푸리에 변환이 실수에 적용될때는 대칭으로 결과가 나오기 때문에, (x [0] == x [x.lenth-1]) 그 반의 데이터(0 to N/2 -1) 그래프를 그렸습니다.
또한 각 인덱스 i에 따른 freqency 계산은 다음과 같습니다.
freq = i * Fs / N;Fs는 sampling frequency이고, N은 FFT 사이즈 입니다. (1024)
N=1024, 샘플링 주기 T=0.004(250Hz)를 가지는 24Hz와 97Hz의 sin파의 합을 fft 변환 하였습니다.
를 fft 하면 진동수 fHz에서 magnitude가 나타나게 됩니다.
그러므로,
인 x는 24Hz와 97Hz의 사인파 합이므로 밑과 같이 fft magnitude그래프가 그려지게 됩니다.
N/2 이후로는 대칭으로 그려지기 때문에, 생략하였습니다.
24와 97은 보기좋게 하기위해 임의로 선택한 숫자 입니다.
밑의 전체 코드를 참조하세요.
//fft
double y[] = new double[N]; //Imaginary Part
for(int i=0;i<N;i++){
y[i] = 0;
} //Imaginary Part는 0으로 채운다.
double x[] = new double[N]; //Real Part
for(int i=0;i<N;i++){
x[i] = Math.sin(2*Math.PI*24*0.004*i) + Math.sin(2*Math.PI*97*0.004*i);
}
FFT doFFT = new FFT(N);
doFFT.fft(x,y); //FFT
double[] mag = new double[N/2]; //magnitude
for(int k=0;k<N/2;k++){
mag[k] = Math.sqrt(Math.pow(x[k],2)+Math.pow(y[k],2));
}
JTransform 라이브러리
Java에서 FFT를 수행하기 위한 유명한 라이브러리 JTransform입니다.
결과는 위와 같은 것을 확인했습니다.
이 라이브러리를 사용하면 1차원의 FFT뿐 아니라 2차원의 FFT(이미지 등)를 할 수 있습니다.
Gradle 추가
dependencies {
...
implementation 'com.github.wendykierp:JTransforms:3.1'
...
}안드로이드에서는 build.gradle(Module:app)에 위와 같이 추가하면 됩니다.
Maven central은 다음과 같습니다.
<dependency>
<groupId>com.github.wendykierp</groupId>
<artifactId>JTransforms</artifactId>
<version>3.1</version>
<classifier>with-dependencies</classifier>
</dependency>
사용 방법
JTransform 라이브러리는 Input과 Output에 실수랑 허수가 순서대로 들어갑니다.
즉, 실수는 배열의 짝수부분, 허수는 배열의 홀수 부분에 들어있다고 생각하면 됩니다.
따라서 사이즈도 그 2배를 가지는 배열을 준비하여야 합니다.
double y[] = new double[N]; //Imaginary Part
for(int i=0;i<N;i++){
y[i] = 0;
}
double x[] = new double[N]; //Real Part
for(int i=0;i<N;i++){
x[i] = Math.sin(2*Math.PI*24*0.004*i) + Math.sin(2*Math.PI*97*0.004*i);
}
double[] a = new double[2*N]; //fft 수행할 배열 사이즈 2N
for(int k=0;k<N;k++){
a[2*k] = x[k]; //Re
a[2*k+1] = y[k]; //Im
}
DoubleFFT_1D fft = new DoubleFFT_1D(N); //1차원의 fft 수행
fft.complexForward(a); //a 배열에 output overwrite
double[] mag = new double[N/2];
for(int k=0;k<N/2;k++){
mag[k] = Math.sqrt(Math.pow(a[2*k],2)+Math.pow(a[2*k+1],2));
}
허수 부분을 0으로 가지는 위와 같은 사인파를 넣었습니다. 보통 전압은 허수부분을 0으로 가지는 실수입니다.
JTransform은 실제 데이터에 complexForward()대신 realForward(), realForwradFull() 메서드를 제공하고 있지만, 이는 complexForward()함수 허수부분을 0으로 처리해 반만 리턴하거나,(대칭이기 때문에) 전부 리턴한것의 차이입니다.
또한 이미지를 FFT 할경우 DoubleFFT_1D 대신 DoubleFFT_2D를 사용하면 됩니다.
FFT class From WikiJava
또 다른 코드를 소개합니다. 사용법은 JTransform 라이브러리랑 비슷합니다.
public class FFTbase {
/**
* The Fast Fourier Transform (generic version, with NO optimizations).
*
* @param inputReal
* an array of length n, the real part
* @param inputImag
* an array of length n, the imaginary part
* @param DIRECT
* TRUE = direct transform, FALSE = inverse transform
* @return a new array of length 2n
*/
public double[] fft(final double[] inputReal, double[] inputImag,
boolean DIRECT) {
// - n is the dimension of the problem
// - nu is its logarithm in base e
int n = inputReal.length;
// If n is a power of 2, then ld is an integer (_without_ decimals)
double ld = Math.log(n) / Math.log(2.0);
// Here I check if n is a power of 2. If exist decimals in ld, I quit
// from the function returning null.
if (((int) ld) - ld != 0) {
System.out.println("The number of elements is not a power of 2.");
return null;
}
// Declaration and initialization of the variables
// ld should be an integer, actually, so I don't lose any information in
// the cast
int nu = (int) ld;
int n2 = n / 2;
int nu1 = nu - 1;
double[] xReal = new double[n];
double[] xImag = new double[n];
double tReal, tImag, p, arg, c, s;
// Here I check if I'm going to do the direct transform or the inverse
// transform.
double constant;
if (DIRECT)
constant = -2 * Math.PI;
else
constant = 2 * Math.PI;
// I don't want to overwrite the input arrays, so here I copy them. This
// choice adds \Theta(2n) to the complexity.
for (int i = 0; i < n; i++) {
xReal[i] = inputReal[i];
xImag[i] = inputImag[i];
}
// First phase - calculation
int k = 0;
for (int l = 1; l <= nu; l++) {
while (k < n) {
for (int i = 1; i <= n2; i++) {
p = bitreverseReference(k >> nu1, nu);
// direct FFT or inverse FFT
arg = constant * p / n;
c = Math.cos(arg);
s = Math.sin(arg);
tReal = xReal[k + n2] * c + xImag[k + n2] * s;
tImag = xImag[k + n2] * c - xReal[k + n2] * s;
xReal[k + n2] = xReal[k] - tReal;
xImag[k + n2] = xImag[k] - tImag;
xReal[k] += tReal;
xImag[k] += tImag;
k++;
}
k += n2;
}
k = 0;
nu1--;
n2 /= 2;
}
// Second phase - recombination
k = 0;
int r;
while (k < n) {
r = bitreverseReference(k, nu);
if (r > k) {
tReal = xReal[k];
tImag = xImag[k];
xReal[k] = xReal[r];
xImag[k] = xImag[r];
xReal[r] = tReal;
xImag[r] = tImag;
}
k++;
}
// Here I have to mix xReal and xImag to have an array (yes, it should
// be possible to do this stuff in the earlier parts of the code, but
// it's here to readibility).
double[] newArray = new double[xReal.length * 2];
double radice = 1 / Math.sqrt(n);
for (int i = 0; i < newArray.length; i += 2) {
int i2 = i / 2;
// I used Stephen Wolfram's Mathematica as a reference so I'm going
// to normalize the output while I'm copying the elements.
newArray[i] = xReal[i2] * radice;
newArray[i + 1] = xImag[i2] * radice;
}
return newArray;
}
/**
* The reference bitreverse function.
*/
private int bitreverseReference(int j, int nu) {
int j2;
int j1 = j;
int k = 0;
for (int i = 1; i <= nu; i++) {
j2 = j1 / 2;
k = 2 * k + j1 - 2 * j2;
j1 = j2;
}
return k;
}
}
사용 코드는 JTransform과 거의 비슷하고, 다만 fft 리턴 값을 2N개의 배열로 받게 됩니다. (점점 설명 짧아지죠..ㅎ)
double x[] = new double[N]; //Real Part
for(int i=0;i<N;i++){
x[i] = Math.sin(2*Math.PI*24*0.004*i) + Math.sin(2*Math.PI*97*0.004*i);
}
double[] y = new double[N]; //Imaginary Part
for(int i=0;i<N;i++) y[i] = 0;
double[] a = new double[2*N];
FFTbase dofft = new FFTbase();
a = dofft.fft(x,y,true);
for(int k=0;k<N/2;k++){
mag[k] = Math.sqrt(Math.pow(a[2*k],2)+Math.pow(a[2*k+1],2));
}
3개의 FFT 코드 비교
3개 각각의 코드를 실제 근전도 신호에 적용해보았습니다.
샘플링 주파수 1500Hz 에 N=16384 (2의 14제곱..) 입니다.
직접 N개의 데이터 엑셀파일을 불러와 각각 FFT 하였습니다.
차례대로 Columbia Univ, JTransform, WikiJava FFT입니다. (소스는 후에 깃헙에 올릴지 생각중..)
실제 데이터에 수행해본 결과, 크게 다른점이 없었습니다.
앞의 2개의 코드는 다른점이 거의 없었고 마지막 코드가 y축의 값이 다르지만 모양은 비슷하였습니다.
다만 dB(Power)를 구햇을 때는 마지막 코드에서 다른점이 많이 보였습니다.
또한 앞의 두개의 코드 그래프를 비교하여 JTransform을 사용하기로 하였습니다. 둘다 0(DC인듯.)에서 큰 값을 보였지만 JTransform은 단순히 0에서의 값을 0으로 처리하면 그래프가 깨끗해지는 양상을 보였습니다.
참고 문서
- FFT library in adroid Sdk
- How do I obtain the frequencies of each value in an FFT?
- JTransform
- JTransform - DoubleFFT_1D
- Columbia University FFT
- FFT From WikiJava
'Android > 통신' 카테고리의 다른 글
| [안드로이드] 아두이노와 안드로이드 Bluetooth 통신하기 (24) | 2020.05.08 |
|---|---|
| [안드로이드 Java] Bandpass Filter(BPF) 구현하기 ― BPF IIR library, FIR 코드 (2) | 2020.03.04 |
| [안드로이드] 프래그먼트(Fragment)에 데이터 전달 하는 법 ― EventBus사용하여 데이터 전달하기 (2) | 2020.02.28 |
| [안드로이드 Java] 하드웨어 모듈과 TCP/IP 소켓 통신 하기 ― 클라이언트 소켓 프로그래밍 (5) | 2020.02.12 |
