16bit FFT

https://gist.github.com/Tomwi/3842231

fix_fft.c

/*   FIX_MPY() - fixed-point multiplication & scaling.   Substitute inline assembly for hardware-specific   optimization suited to a particluar DSP processor.   Scaling ensures that result remains 16-bit. */ inline short FIX_MPY(short a,short b) {     /* shift right one less bit (i.e. 15-1) */      int c = ((int)a * (int)b) >> 14;      /* last bit shifted out = rounding-bit */      b = c & 0x01;      /* last shift + rounding bit */      a = (c >> 1) + b;      return a; }

 
이 코드는 고정 소수점(fixed-point) 형식을 사용하여 두 16비트 정수(short)의 곱셈을 수행하는 함수인 FIX_MPY()를 정의합니다. 이 함수는 결과가 여전히 16비트로 유지되도록 스케일링합니다.

주요 요소는 다음과 같습니다:

1. FIX_MPY 함수
두 개의 16비트 정수를 입력으로 받아서 곱셈을 수행하고, 결과를 16비트로 반환하는 함수입니다.

2. 스케일링
스케일링은 고정 소수점 연산에서 중요한 요소입니다. 여기서는 결과가 16비트로 유지되도록 스케일링됩니다.

3. 비트 시프트 연산
곱셈 결과를 계산하기 위해 두 정수를 int(32비트) 형식으로 캐스팅한 후, 오른쪽으로 14비트 시프트합니다. 이렇게 함으로써 결과의 정수 부분은 16비트에서 15비트로 줄어들게 되며, 그 결과는 스케일링된 값이 됩니다.

4. 라운딩
곱셈 결과를 계산한 후에는 라운딩이 수행됩니다. 이는 마지막 비트가 반올림 비트가 되도록 하기 위한 것입니다. 즉, 최하위 비트를 0 또는 1로 설정하여 결과를 반올림합니다.

5. 결과 반환
최종 결과는 16비트 정수로 반환됩니다.

이 함수는 다양한 DSP(Digital Signal Processor) 프로세서에 맞게 인라인 어셈블리로 대체되어야 하며, 특정 하드웨어에 최적화될 수 있습니다. 이 코드는 플랫폼에 특정하지 않고 범용적으로 사용할 수 있는 것을 목표로 하고 있습니다.
 
 
 
 
 

/*   fix_fft() - perform forward/inverse fast Fourier transform.   fr[n],fi[n] are real and imaginary arrays,both INPUT AND   RESULT (in-place FFT),with 0 <= n < 2**m; set inverse to   0 for forward transform (FFT),or 1 for iFFT. */ int fix_fft(short fr[],short fi[],short m,short inverse) {     int mr,nn,i,j,l,k,istep,n,scale,shift;      short qr,qi,tr,ti,wr,wi;      n = 1 << m;      /* max FFT size = N_WAVE */      if (n > N_WAVE)           return -1;      mr = 0;      nn = n - 1;      scale = 0;      /* decimation in time - re-order data */      for (m=1; m<=nn; ++m) {           l = n;           do {                l >>= 1;           } while (mr+l > nn);           mr = (mr & (l-1)) + l;           if (mr <= m)                continue;           tr = fr[m];           fr[m] = fr[mr];           fr[mr] = tr;           ti = fi[m];           fi[m] = fi[mr];           fi[mr] = ti;      }      l = 1;      k = LOG2_N_WAVE-1;      while (l < n) {           if (inverse) {                /* variable scaling,depending upon data */                shift = 0;                for (i=0; i<n; ++i) {                     j = fr[i];                     if (j < 0)                          j = -j;                     m = fi[i];                     if (m < 0)                          m = -m;                     if (j > 16383 || m > 16383) {                          shift = 1;                          break;                     }                }                if (shift)                     ++scale;                } else {                     /*                       fixed scaling,for proper normalization --                       there will be log2(n) passes,so this results                       in an overall factor of 1/n,distributed to                       maximize arithmetic accuracy.                     */                     shift = 1;                }                /*                  it may not be obvious,but the shift will be                  performed on each data point exactly once,                  during this pass.                */                istep = l << 1;                for (m=0; m<l; ++m) {                     j = m << k;                     /*  0 <= j < N_WAVE/2 */                     wr =  Sinewave[j+N_WAVE/4];                     wi = -Sinewave[j];                     if (inverse)                          wi = -wi;                     if (shift) {                          wr >>= 1;                          wi >>= 1;                     }                for (i=m; i<n; i+=istep) {                     j = i + l;                     tr = FIX_MPY(wr,fr[j]) - FIX_MPY(wi,fi[j]);                     ti = FIX_MPY(wr,fi[j]) + FIX_MPY(wi,fr[j]);                     qr = fr[i];                     qi = fi[i];                     if (shift) {                          qr >>= 1;                          qi >>= 1;                     }                     fr[j] = qr - tr;                     fi[j] = qi - ti;                     fr[i] = qr + tr;                     fi[i] = qi + ti;                }           }           --k;           l = istep;      }      return scale; }

이 코드는 고정 소수점(fixed-point) FFT(Fast Fourier Transform) 알고리즘인 fix_fft() 함수를 구현한 것입니다. FFT는 주파수 영역으로 신호를 변환하는 방법으로, 신호 처리나 데이터 압축 등 다양한 응용 분야에서 사용됩니다. fix_fft() 함수는 입력으로 주어진 실수 배열인 fr[]과 허수 배열인 fi[]에 대해 FFT 또는 IFFT(Inverse FFT)를 수행합니다.

주요 내용을 간략히 설명하겠습니다:

1. 함수 시그니처
함수는 입력으로 실수 배열 fr[], 허수 배열 fi[], FFT의 크기를 나타내는 m, 그리고 forward transform (FFT)를 할 것인지 inverse transform (iFFT)를 할 것인지 나타내는 inverse를 받습니다.
예) 데이터 개수는 2의 배수여야하고, 1024개일때 m=9입니다. 2^10=1024

2. 주요 변수 및 초기화
함수 내에서 사용되는 변수들을 선언하고 초기화합니다.

3. 데이터 재배열
주어진 데이터를 재배열합니다. 이는 FFT 알고리즘을 위한 입력 데이터 구조를 만들기 위한 것입니다.

4. FFT 또는 IFFT 수행
FFT 또는 IFFT를 수행합니다. 이는 주어진 데이터에 대해 주파수 영역으로의 변환을 의미합니다. 이 과정은 푸리에 변환(Fourier Transform)을 효율적으로 수행하기 위한 방법인 Cooley-Tukey 알고리즘을 기반으로 합니다.

5. 스케일링
결과를 스케일링합니다. 스케일링은 주로 IFFT를 수행할 때 사용되며, 정규화를 위해 데이터에 적절한 스케일링을 적용합니다.

6. 결과 반환
변환된 데이터를 반환하고, 필요에 따라 스케일링된 횟수를 반환합니다.

이 코드는 FFT 알고리즘을 효율적으로 구현하기 위해 고정 소수점 연산을 사용하며, 주어진 크기의 FFT에 대한 연산을 수행합니다.
 
 

/*   fix_fftr() - forward/inverse FFT on array of real numbers.   Real FFT/iFFT using half-size complex FFT by distributing   even/odd samples into real/imaginary arrays respectively.   In order to save data space (i.e. to avoid two arrays,one   for real,one for imaginary samples),we proceed in the   following two steps: a) samples are rearranged in the real   array so that all even samples are in places 0-(N/2-1) and   all imaginary samples in places (N/2)-(N-1),and b) fix_fft   is called with fr and fi pointing to index 0 and index N/2   respectively in the original array. The above guarantees   that fix_fft "sees" consecutive real samples as alternating   real and imaginary samples in the complex array. */ int fix_fftr(short f[],  int m,//m=9(1024샘플)  int inverse) {     int i,N = 1<<(m-1),scale = 0;      short tt,*fr=f,*fi=&f[N];      if (inverse)           scale = fix_fft(fi,fr,m-1,inverse);      for (i=1; i<N; i+=2) {           tt = f[N+i-1];           f[N+i-1] = f[i];           f[i] = tt;      }      if (! inverse)           scale = fix_fft(fi,fr,m-1,inverse);      return scale; }

위의 코드는 실수 배열에 대한 FFT(Fast Fourier Transform) 및 IFFT(Inverse Fast Fourier Transform)를 수행하는 함수인 fix_fftr()을 구현한 것입니다. 
이 함수는 주어진 실수 배열에서 FFT 또는 IFFT를 수행하며, 배열의 크기가 2의 거듭제곱이어야 합니다.

1. 함수 시그니처
함수는 입력으로 실수 배열 f[], 배열의 크기를 나타내는 m, 그리고 forward transform (FFT)를 할 것인지 inverse transform (iFFT)를 할 것인지 나타내는 inverse를 받습니다.

2. 주요 변수 및 초기화
함수 내에서 사용되는 변수들을 선언하고 초기화합니다. N은 FFT 또는 IFFT의 크기를 나타내는 변수입니다.

3. 배열 재배열
주어진 데이터를 재배열하여 모든 짝수 인덱스는 실수 부분에, 홀수 인덱스는 허수 부분에 배치합니다. 이렇게 함으로써 실수 배열을 복소수 배열로 변환하여 FFT 또는 IFFT를 적용할 수 있습니다.

4. FFT 또는 IFFT 수행
FFT 또는 IFFT를 수행합니다.
배열의 실수 부분(fr)과 허수 부분(fi)을 바꿔서 fix_fft() 함수를 호출합니다.

5. 스케일링
결과를 스케일링합니다. 스케일링은 FFT 또는 IFFT를 수행할 때 사용됩니다.
결과의 크기를 조정하여 정규화를 수행합니다.

6. 결과 반환
변환된 데이터를 반환하고, 필요에 따라 스케일링된 횟수를 반환합니다.

이 함수는 실수 배열에 대한 FFT 또는 IFFT를 수행하는데, 입력 데이터를 복소수 형태로 재배열하여 FFT를 적용합니다. 
이렇게 함으로써 복소수 배열에 대한 FFT를 사용하여 실수 배열에 대한 FFT를 수행할 수 있습니다.
 
 

정규화

FFT 결과를 정규화하는 것은 변환된 주파수 영역 데이터를 원래의 시간 영역 데이터와 동등한 단위로 표현하는 과정입니다. 주로 두 가지 방법으로 정규화를 수행합니다.

1. 에너지 정규화(Energy normalization)
   이 방법은 FFT 결과의 크기를 원래 입력 신호의 에너지와 일치시킵니다. 주파수 영역에서 각 항목의 크기를 입력 신호의 에너지에 비례하도록 조정합니다. 이를 위해 FFT 결과를 절대값으로 취하고, 스케일링 상수를 사용하여 크기를 조정합니다.

2. 피크 정규화(Peak normalization)
   이 방법은 FFT 결과에서 가장 큰 값을 1로 정규화합니다. 가장 큰 값에 해당하는 주파수의 크기를 1로 만들고, 다른 주파수 값들을 이에 상대적으로 조정합니다. 이 방법은 각 주파수의 상대적인 크기를 비교할 때 유용합니다.

정규화를 통해 FFT 결과를 일반화하고, 사용자가 결과를 해석하거나 다른 처리를 수행하기 쉽도록 만듭니다. 주파수 영역 데이터의 정확한 크기나 상대적인 크기를 비교하는 등의 작업을 수행할 때 유용합니다.
 

스케일링을 하는 이유

FFT(Fast Fourier Transform)를 수행하는 과정에서 스케일링을 하는 이유는 주로 다음과 같습니다:

1. 정규화
FFT를 수행하면 주파수 영역으로의 변환을 얻게 됩니다. 이 변환된 결과는 원래의 입력 신호와 크기 및 단위가 다르기 때문에 정규화가 필요합니다. 스케일링을 통해 결과를 정규화하여 사용자가 이해하기 쉽고 특정 응용 분야에 맞게 조정할 수 있습니다.

2. 동적 범위 조정
FFT 결과의 동적 범위가 너무 넓을 수 있습니다. 이는 결과가 너무 커서 계산 또는 표시하기가 어려울 수 있습니다. 스케일링을 통해 결과의 동적 범위를 적절하게 조정하여 계산 및 표시를 용이하게 할 수 있습니다.

3. 정밀도와 오차 제어
고정 소수점 연산을 사용하는 경우, 결과의 정밀도와 오차를 관리하는 것이 중요합니다. 스케일링을 통해 결과의 정밀도를 유지하고 오차를 최소화할 수 있습니다.

4. 변환된 데이터의 해석
FFT 결과를 해석할 때 스케일링은 결과를 해석하기 쉽게 만들어 줍니다. 예를 들어, 주파수 영역에서의 특정 구조나 패턴을 쉽게 식별할 수 있도록 도와줍니다.

따라서 FFT 결과를 실제로 활용하기 위해서는 결과를 적절하게 스케일링하여 정규화하고, 계산 또는 표시하는 데 필요한 형식으로 조정하는 것이 중요합니다.
 
 

오디오 16비트 signed 데이터를 넣는 방법

위의 코드에서 오디오 16비트 signed 데이터를 사용할 때, 오디오 데이터를 그대로 넣으면 됩니다.
주어진 함수는 입력으로 실수 배열 `fr[]`과 허수 배열 `fi[]`를 받는데, 여기서 `fr[]` 배열에는 오디오의 실수 부분이, `fi[]` 배열에는 오디오의 허수 부분이 저장될 것입니다.

오디오 데이터를 FFT에 사용할 때 주의해야 할 몇 가지 사항이 있습니다.

1. 입력 데이터의 형식
오디오 데이터는 16비트 signed 형식으로 표현됩니다. 이는 정수로 표현되며, 주로 -32768에서 32767까지의 범위를 가집니다.

2. 허수 부분의 초기화
대부분의 경우 오디오 데이터는 실수 값만을 포함하므로, 허수 부분(`fi[]` 배열)은 일반적으로 0으로 초기화합니다. 오디오 데이터를 FFT에 사용할 때 허수 부분은 대개 사용되지 않습니다.

따라서 오디오 데이터의 실수 부분을 `fr[]` 배열에 넣고, `fi[]` 배열은 0으로 초기화하여 함수에 전달하면 됩니다. 이렇게 하면 FFT가 정상적으로 수행될 것입니다.
 
 
 

오디오파일 FFT

44.1kHz singed 16bit mono wav파일데이터를 한번에 1024개씩 FFT하면,
1024개의 실수 부분(fr)과 허수 부분(fi)값이 나오며, 
스펙트럼 아날라이져로 그리기위해서는 아래의 값을 그려줍니다.
sqrt(실수 부분(fr)*실수 부분(fr) + 허수 부분(fi)* 허수 부분(fi))



인터넷에 올라온 fft소스코드중에 ifft까지 잘되는것 찾기가 쉽지않음.