FFT, Fast Fourier transform. soft for “C” compiler
Довелось в програмі для мікроконтролера обчислювати спектр оцифрованого сигналу. Досить швидко я переконався, що теорії в Інтернеті можна знайти багато, але хочеться мати готовий програмний зразок, взяв – і використав.
Щоб хтось інший не повторював пошук і оптимізацію програми для обчислення спектра, всім бажаючим надаю готову, перевірену підпрограму, я використовував її для одного з мікроконтролерів, вона написана на “C”.
Трошки теорії.
Для обчислення спектра беремо виборку з Nmas оцифрованих точок якогось сигналу, Nmas – це число, яке дорівнює “2 в якомусь ступені”. Наприклад, Nmas = 512, 1024, 2048, і так далі.
Кількість точок спектра, які ми отримаємо, в два рази менша.
Кожна точка спектра представлена двома числами, це дійсна і уявна частина.
Якщо сигнал дискретизований з частотою F, то одержуємо спектр в частотній полосі від 0 до F/2.
Якщо бездумно порахувати спектр, він буде поганенький, бо наша виборка не нескінченна. Щоб покращити результат, використовується нормуюча функція. Якою вона має бути, спеціалісти вже порахували до нас. Просто використовуємо.
—————————–
отже, підпрограма
—————————–
// Необхідні величини:
#define Nmas 1024 // кількість точок виборки для вирахування спектра (наприклад, 1024)
#define RozmREZ 2048 // масив для результатів, попарно дійсна і уявна частина. RozmREZ=Nmas*2
#define ForNorm 513 // в два рази менше від Nmas + 1 (Nmas/2+1)
#define PorMas 10 // ступінь числа 2 для одержання Nmas (двійковий логарифм від Nmas)
// наприклад:
// PorMas=1 Nmas=2
// PorMas=2 Nmas=4
// PorMas=3 Nmas=8
// PorMas=4 Nmas=16
// PorMas=5 Nmas=32
// PorMas=6 Nmas=64
// PorMas=7 Nmas=128
// PorMas=8 Nmas=256
// PorMas=9 Nmas=512
// PorMas=10 Nmas=1024
// і так далі
// використовуються масиви:
float AVAL[Nmas]; // вхідні дані
float REZ[Nmas+Nmas]; // результат
void FFT_Prog(void)
{
//
// Обчислення спектра, швидке перетворення Фур’є
// FFT, Fast Fourier transform
//
// AVAL – масив вхідних данних, Nmas – довжина масива, це 2 в вибраному вамі ступені, наприклад, 2^10=1024 .
//
// REZ – масив одержаних значень, попарно дійсна і уявна частина,
// при цьому нам корисна в масиві лише перша половина пар,
// бо інша половина дзеркально подібна на першу.
// REZ[0], REZ[1] – перша точка спектра, дійсна і уявна частина.
// REZ[2], REZ[3] – друга точка спектра, дійсна і уявна частина, і так далі.
// Наприклад, при Nmas=1024 у нас є 512 точок спектра, кожна точка
// представлена двома числами, дійсною і уявною частиною.
// Розмір масиву REZ в два рази більший, ніж масиву AVAL
// (але друга половина масиву буде дзеркальною до першої, і нам не потрібна)
//
// Використовується одноразовий попередній прорахунок
// нормуючої функції, щоб кожен раз при звертанні до функції не обчислювати.
// Підпрограма Фур’є робить обчислення лише (log2(Nmax))*2 синусів,
// тобто, наприклад, при Nmax=1024 робиться 10*2=20 обчислень синусів
// (навіть їх обчислення робиться один перший раз).
// Попереднє одноразове обчислення збільшує наступну швидкість алгоритму.
//
// стартові дані беремо з масиву AVAL
// результат обчислень записано в масиві REZ
//
static int priznorm=0, indnor, napnor, Nnor;
//
// В масиві NORM – один раз розрахована нормуюча функція
// в вигляді масива потрібних значень
// Признак “static” в підпрограмі є обов’язковим
//
static float NORM[ForNorm]; // розмір масива = Nmas/2 + 1
//
// один раз розраховані потрібні синуси MasWp і MasWtmp
// (щоб кожен раз не рахувати)
//
static float MasWpi[PorMas]; // розмір масива = PorMas
static float MasWtmp[PorMas]; // розмір масива = PorMas
//
// лічильники і тимчасові змінні
//
int i, j, N, M, Mmax, Istp, CounSi;
float Tmpr, Tmpi, Wtmp, Theta;
float Wpr, Wpi, Wr, Wi;
//
// початок обчислень
//
Nnor = Nmas/2;
N=2;
if(priznorm==0){
//
// попередній одноразовий прорахунок необхідних синусів
//
for (i = 0; i < PorMas; i+=1) {
Theta = -TwoPi / N; MasWpi[i] = sin(Theta);
MasWtmp[i] = sin(Theta / 2); N=N*2; };
//
// попередній одноразовий прорахунок нормуючої функції
//
for (i = 0; i <= Nnor; i+=1) {
NORM[i] = 0.54-0.46*cos(i*TwoPi/Nmas);}; priznorm=1;};
//
// Це лише одна з нормуючих функцій, з коефіцієнтами 0.54 і 0.46,
// вона вважається досить вдалою
//
// Спочатку формуємо початковий масив пронормованих даних
//
N = Nmas * 2;
indnor=0;napnor=0;
for (i = 0; i < N; i+=2) {
REZ[i] = 0.0;
REZ[i+1] = AVAL[i/2]*NORM[indnor]; // нормування віконною функцією
if(napnor==0)
{indnor++;if(indnor==Nnor)napnor=1;}
else{indnor–;if(indnor==0)napnor=0;};
};
//
// Зформували масив пронормованих даних.
// Цей масив прораховується один раз.
// Далі – саме обчислення швидкого перетворення Фур’є
//
i = 1; j = 1;
while (i < N) {
if (j > i) {
Tmpr = REZ[i]; REZ[i] = REZ[j]; REZ[j] = Tmpr;
Tmpr = REZ[i+1]; REZ[i+1] = REZ[j+1]; REZ[j+1] = Tmpr;
};
i = i + 2; M = Nmas;
while ((M >= 2) && (j > M)) {
j = j – M; M = M >> 1;
};
j = j + M;
};
Mmax = 2; CounSi=0;
while (N > Mmax) {
//
// Theta = -TwoPi / Mmax; // могло бути звичайне обчислення
// Wpi = sin(Theta); // синусів
// Wtmp = sin(Theta / 2);
Wpi = MasWpi[CounSi]; // використовуємо готове обчислене значення
Wtmp = MasWtmp[CounSi];
//
Wpr = Wtmp * Wtmp * 2;
Istp = Mmax * 2; Wr = 1; Wi = 0; M = 1;
while (M < Mmax) {
i = M; M = M + 2; Tmpr = Wr; Tmpi = Wi;
Wr = Wr – Tmpr * Wpr – Tmpi * Wpi;
Wi = Wi + Tmpr * Wpi – Tmpi * Wpr;
while (i < N) {
j = i + Mmax;
Tmpr = Wr * REZ[j] – Wi * REZ[j-1];
Tmpi = Wi * REZ[j] + Wr * REZ[j-1];
REZ[j] = REZ[i] – Tmpr; REZ[j-1] = REZ[i-1] – Tmpi;
REZ[i] = REZ[i] + Tmpr; REZ[i-1] = REZ[i-1] + Tmpi;
i = i + Istp;
};
};
Mmax = Istp; CounSi++;
};
}
Весь список статей нашого сайту, корисних для програмістів, знаходиться ось тут