Openframeworksでビジュアライズ

VJで使ってるFFT演算のサンプルです。
Openframeworksにはマイク入力の音量を取得するためのツールがあんま無いので、以下のフォーラムのコードを拝借してやっています。
Forum:examples ~ audio input FFT example





まずヘッダの中身

[testApp.h]

#ifndef _TEST_APP
#define _TEST_APP

#define BUFFER_SIZE 1024
#define NUM_WINDOWS 80

#include "ofMain.h"
#include "sound_input.h"

class testApp : public ofBaseApp{

public:
void setup();
void update();
void draw();

void keyPressed  (int key);
void keyReleased(int key);
void mouseMoved(int x, int y );
void mouseDragged(int x, int y, int button);
void mousePressed(int x, int y, int button);
void mouseReleased(int x, int y, int button);
void windowResized(int w, int h);

void sound_init();
void sound_get();
void audioReceived  (float * input, int bufferSize, int nChannels); 
float * left;
float * right;
float gain;
float level;
float timeline_array[257][2];
int  bufferCounter;
fft  myfft;
float phase[BUFFER_SIZE];
float power[BUFFER_SIZE];
float freq[NUM_WINDOWS][BUFFER_SIZE/2];
float freq_phase[NUM_WINDOWS][BUFFER_SIZE/2];
float magnitude[BUFFER_SIZE];

};

#endif

#defineと、段落が下がってる関数・変数が追加されています。それとこの後追加するsound_input.hもインクルードする必要があります。

音取るだけなんですが、計算に結構な数のパラメータが必要です。
そして別途sound_input.h及びsound_input.cppをプロジェクトに追加します。

[sound_input.h]

#ifndef _FFT
#define _FFT

#ifndef M_PI
#define M_PI  3.14159265358979323846  /* pi */
#endif

class fft {

public:

fft();
~fft(); 

/* Calculate the power spectrum */
void powerSpectrum(int start, int half, float *data, int windowSize,float *magnitude,float *phase, float *power, float *avg_power);
/* ... the inverse */
void inversePowerSpectrum(int start, int half, int windowSize, float *finalOut,float *magnitude,float *phase);  

};


#endif

[sound_input.cpp]

#include "sound_input.h"

#include 
#include 
#include 

int **gFFTBitTable = NULL;
const int MaxFastBits = 16;

int IsPowerOfTwo(int x)
{
if (x < 2)
return false;

if (x & (x - 1))             /* Thanks to 'byang' for this cute trick! */
return false;

return true;
}

int NumberOfBitsNeeded(int PowerOfTwo)
{
int i;

if (PowerOfTwo < 2) {
fprintf(stderr, "Error: FFT called with size %d\n", PowerOfTwo);
exit(1);
}

for (i = 0;; i++)
if (PowerOfTwo & (1 << i))
return i;
}

int ReverseBits(int index, int NumBits)
{
int i, rev;

for (i = rev = 0; i < NumBits; i++) {
rev = (rev << 1) | (index & 1);
index >>= 1;
}

return rev;
}

void InitFFT()
{
gFFTBitTable = new int *[MaxFastBits];

int len = 2;
for (int b = 1; b <= MaxFastBits; b++) {

gFFTBitTable[b - 1] = new int[len];

for (int i = 0; i < len; i++)
gFFTBitTable[b - 1][i] = ReverseBits(i, b);

len <<= 1;
}
}

inline int FastReverseBits(int i, int NumBits)
{
if (NumBits <= MaxFastBits)
return gFFTBitTable[NumBits - 1][i];
else
return ReverseBits(i, NumBits);
}

/*
* Complex Fast Fourier Transform
*/

void FFT(int NumSamples,
bool InverseTransform,
float *RealIn, float *ImagIn, float *RealOut, float *ImagOut)
{
int NumBits;                 /* Number of bits needed to store indices */
int i, j, k, n;
int BlockSize, BlockEnd;

double angle_numerator = 2.0 * M_PI;
float tr, ti;                /* temp real, temp imaginary */

if (!IsPowerOfTwo(NumSamples)) {
fprintf(stderr, "%d is not a power of two\n", NumSamples);
exit(1);
}

if (!gFFTBitTable)
InitFFT();

if (InverseTransform)
angle_numerator = -angle_numerator;

NumBits = NumberOfBitsNeeded(NumSamples);

/*
**   Do simultaneous data copy and bit-reversal ordering into outputs...
*/

for (i = 0; i < NumSamples; i++) {
j = FastReverseBits(i, NumBits);
RealOut[j] = RealIn[i];
ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
}

/*
**   Do the FFT itself...
*/

BlockEnd = 1;
for (BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {

double delta_angle = angle_numerator / (double) BlockSize;

float sm2 = sin(-2 * delta_angle);
float sm1 = sin(-delta_angle);
float cm2 = cos(-2 * delta_angle);
float cm1 = cos(-delta_angle);
float w = 2 * cm1;
float ar0, ar1, ar2, ai0, ai1, ai2;

for (i = 0; i < NumSamples; i += BlockSize) {
ar2 = cm2;
ar1 = cm1;

ai2 = sm2;
ai1 = sm1;

for (j = i, n = 0; n < BlockEnd; j++, n++) {
ar0 = w * ar1 - ar2;
ar2 = ar1;
ar1 = ar0;

ai0 = w * ai1 - ai2;
ai2 = ai1;
ai1 = ai0;

k = j + BlockEnd;
tr = ar0 * RealOut[k] - ai0 * ImagOut[k];
ti = ar0 * ImagOut[k] + ai0 * RealOut[k];

RealOut[k] = RealOut[j] - tr;
ImagOut[k] = ImagOut[j] - ti;

RealOut[j] += tr;
ImagOut[j] += ti;
}
}

BlockEnd = BlockSize;
}

/*
**   Need to normalize if inverse transform...
*/

if (InverseTransform) {
float denom = (float) NumSamples;

for (i = 0; i < NumSamples; i++) {
RealOut[i] /= denom;
ImagOut[i] /= denom;
}
}
}

/*
* Real Fast Fourier Transform
*
* This function was based on the code in Numerical Recipes in C.
* In Num. Rec., the inner loop is based on a single 1-based array
* of interleaved real and imaginary numbers.  Because we have two
* separate zero-based arrays, our indices are quite different.
* Here is the correspondence between Num. Rec. indices and our indices:
*
* i1  <->  real[i]
* i2  <->  imag[i]
* i3  <->  real[n/2-i]
* i4  <->  imag[n/2-i]
*/

void RealFFT(int NumSamples, float *RealIn, float *RealOut, float *ImagOut)
{
int Half = NumSamples / 2;
int i;

float theta = M_PI / Half;

float *tmpReal = new float[Half];
float *tmpImag = new float[Half];

for (i = 0; i < Half; i++) {
tmpReal[i] = RealIn[2 * i];
tmpImag[i] = RealIn[2 * i + 1];
}

FFT(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);

float wtemp = float (sin(0.5 * theta));

float wpr = -2.0 * wtemp * wtemp;
float wpi = float (sin(theta));
float wr = 1.0 + wpr;
float wi = wpi;

int i3;

float h1r, h1i, h2r, h2i;

for (i = 1; i < Half / 2; i++) {

i3 = Half - i;

h1r = 0.5 * (RealOut[i] + RealOut[i3]);
h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
h2i = -0.5 * (RealOut[i] - RealOut[i3]);

RealOut[i] = h1r + wr * h2r - wi * h2i;
ImagOut[i] = h1i + wr * h2i + wi * h2r;
RealOut[i3] = h1r - wr * h2r + wi * h2i;
ImagOut[i3] = -h1i + wr * h2i + wi * h2r;

wr = (wtemp = wr) * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}

RealOut[0] = (h1r = RealOut[0]) + ImagOut[0];
ImagOut[0] = h1r - ImagOut[0];

delete[]tmpReal;
delete[]tmpImag;
}

/*
* PowerSpectrum
*
* This function computes the same as RealFFT, above, but
* adds the squares of the real and imaginary part of each
* coefficient, extracting the power and throwing away the
* phase.
*
* For speed, it does not call RealFFT, but duplicates some
* of its code.
*/

void PowerSpectrum(int NumSamples, float *In, float *Out)
{
int Half = NumSamples / 2;
int i;

float theta = M_PI / Half;

float *tmpReal = new float[Half];
float *tmpImag = new float[Half];
float *RealOut = new float[Half];
float *ImagOut = new float[Half];

for (i = 0; i < Half; i++) {
tmpReal[i] = In[2 * i];
tmpImag[i] = In[2 * i + 1];
}

FFT(Half, 0, tmpReal, tmpImag, RealOut, ImagOut);

float wtemp = float (sin(0.5 * theta));

float wpr = -2.0 * wtemp * wtemp;
float wpi = float (sin(theta));
float wr = 1.0 + wpr;
float wi = wpi;

int i3;

float h1r, h1i, h2r, h2i, rt, it;
//float total=0;

for (i = 1; i < Half / 2; i++) {

i3 = Half - i;

h1r = 0.5 * (RealOut[i] + RealOut[i3]);
h1i = 0.5 * (ImagOut[i] - ImagOut[i3]);
h2r = 0.5 * (ImagOut[i] + ImagOut[i3]);
h2i = -0.5 * (RealOut[i] - RealOut[i3]);

rt = h1r + wr * h2r - wi * h2i; //printf("Realout%i = %f",i,rt);total+=fabs(rt);
it = h1i + wr * h2i + wi * h2r; // printf("  Imageout%i = %f\n",i,it);

Out[i] = rt * rt + it * it;

rt = h1r - wr * h2r + wi * h2i;
it = -h1i + wr * h2i + wi * h2r;

Out[i3] = rt * rt + it * it;

wr = (wtemp = wr) * wpr - wi * wpi + wr;
wi = wi * wpr + wtemp * wpi + wi;
}
//printf("total = %f\n",total);
rt = (h1r = RealOut[0]) + ImagOut[0];
it = h1r - ImagOut[0];
Out[0] = rt * rt + it * it;

rt = RealOut[Half / 2];
it = ImagOut[Half / 2];
Out[Half / 2] = rt * rt + it * it;

delete[]tmpReal;
delete[]tmpImag;
delete[]RealOut;
delete[]ImagOut;
}

/*
* Windowing Functions
*/

int NumWindowFuncs()
{
return 4;
}

char *WindowFuncName(int whichFunction)
{
switch (whichFunction) {
default:
case 0:
return "Rectangular";
case 1:
return "Bartlett";
case 2:
return "Hamming";
case 3:
return "Hanning";
}
}

void WindowFunc(int whichFunction, int NumSamples, float *in)
{
int i;

if (whichFunction == 1) {
// Bartlett (triangular) window
for (i = 0; i < NumSamples / 2; i++) {
in[i] *= (i / (float) (NumSamples / 2));
in[i + (NumSamples / 2)] *=
(1.0 - (i / (float) (NumSamples / 2)));
}
}

if (whichFunction == 2) {
// Hamming
for (i = 0; i < NumSamples; i++)
in[i] *= 0.54 - 0.46 * cos(2 * M_PI * i / (NumSamples - 1));
}

if (whichFunction == 3) {
// Hanning
for (i = 0; i < NumSamples; i++)
in[i] *= 0.50 - 0.50 * cos(2 * M_PI * i / (NumSamples - 1));
}
}

/* constructor */
fft::fft() {


}

/* destructor */
fft::~fft() {


}

/* Calculate the power spectrum */
void fft::powerSpectrum(int start, int half, float *data, int windowSize,float *magnitude,float *phase, float *power, float *avg_power) {
int i;
int windowFunc = 3;
float total_power = 0.0f;

/* processing variables*/
float *in_real = new float[windowSize];
float *in_img = new float[windowSize];
float *out_real = new float[windowSize];
float *out_img = new float[windowSize];

for (i = 0; i < windowSize; i++) {
in_real[i] = data[start + i];
}

WindowFunc(windowFunc, windowSize, in_real);
RealFFT(windowSize, in_real, out_real, out_img);

for (i = 0; i < half; i++) {
/* compute power */
power[i] = out_real[i]*out_real[i] + out_img[i]*out_img[i];
total_power += power[i];
/* compute magnitude and phase */
magnitude[i] = 2.0*sqrt(power[i]);
phase[i] = atan2(out_img[i],out_real[i]);
}
/* calculate average power */
*(avg_power) = total_power / (float) half;

delete[]in_real;
delete[]in_img;   
delete[]out_real;
delete[]out_img;
}

void fft::inversePowerSpectrum(int start, int half, int windowSize, float *finalOut,float *magnitude,float *phase) {
int i;
int windowFunc = 3;

/* processing variables*/
float *in_real = new float[windowSize];
float *in_img = new float[windowSize];
float *out_real = new float[windowSize];
float *out_img = new float[windowSize];

/* get real and imag part */
for (i = 0; i < half; i++) { 
in_real[i] = magnitude[i]*cos(phase[i]);
in_img[i]  = magnitude[i]*sin(phase[i]);
}

/* zero negative frequencies */
for (i = half; i < windowSize; i++) { 
in_real[i] = 0.0;
in_img[i] = 0.0;
}

FFT(windowSize, 1, in_real, in_img, out_real, out_img); // second parameter indicates inverse transform
WindowFunc(windowFunc, windowSize, out_real);

for (i = 0; i < windowSize; i++) {
finalOut[start + i] += out_real[i];
}

delete[]in_real;
delete[]in_img;   
delete[]out_real;
delete[]out_img;
}

ここについては完全に呪文です。理解できたらサンプルなんて使わなくても自分で組めます。 そして、testAppクラスに、ヘッダで宣言した以下の三つの関数を追加します。 [testApp.cpp]

void testApp::sound_init(){
// 0 output channels, 
// 2 input channels
// 44100 samples per second
// BUFFER_SIZE samples per buffer
// 4 num buffers (latency)

ofSoundStreamSetup(0,2,this, 44100,BUFFER_SIZE, 4); 
left = new float[BUFFER_SIZE];
right = new float[BUFFER_SIZE];

for (int i = 0; i < NUM_WINDOWS; i++)
{
for (int j = 0; j < BUFFER_SIZE/2; j++)
{
freq[i][j] = 0; 
}
}
}


void testApp::sound_get(){
static int index=0;
float avg_power = 0.0f; 

if(index < 80)
index += 1;
else
index = 0;

/* do the FFT */
myfft.powerSpectrum(0,(int)BUFFER_SIZE/2, left,BUFFER_SIZE,&magnitude[0],&phase[0],&power[0],&avg_power);

/* start from 1 because mag[0] = DC component */
/* and discard the upper half of the buffer */
level = 0;
for(int j=1; j < BUFFER_SIZE/2; j++) {
freq[index][j] = magnitude[j]*gain*3;
level += magnitude[j];
magnitude[j] *= gain*2;
}
level = level*gain;
}


void testApp::audioReceived(float * input, int bufferSize, int nChannels){ 
// samples are "interleaved"
for (int i = 0; i < bufferSize; i++){
left[i] = input[i*2];
right[i] = input[i*2+1];
}
bufferCounter++;
}

これで音を取る準備が完了です。 testAppのSetupとUpdateに以下の文を付け加えます

void testApp::setup(){
gain = 1;
sound_init();
ofSetBackground(0,0,0);
}

void testApp::update(){
sound_get();
}

そして、draw関数内にこんな感じで

void testApp::draw(){ 
//音の大きさに合わせて円を描画
ofNoFill();
ofSetColor(0xFFFFFF);
ofCircle(ofGetWidth()/2, ofGetHeight()/2, level*0.1f);

//周波数スペクトラムの描画
ofSetColor(0xFF0000);
for (int i = 0;i < BUFFER_SIZE;i++){
ofLine(i,400, i, 400-magnitude[i]);
}
}

levelが音量、magnitude[0〜BUFFER_SIZE]の配列がFFTによる各周波数の大きさを表しています。
sound_get()を追うと分かりますが、gainの値で入力の大きさも変えることができます。

ただ元の音量が低すぎると精度が低いので、システム環境設定のサウンドで上げられる所まで上げた方が良いです。
ちなみに入力ソースは、起動した瞬間サウンド設定で選択されてる入力機器が選ばれるようになってるようです。

この状態でlevelとmagnitudeを色んな描画の引数に当てはめると、面白い動きとかします。お試しアレ。