Быстрое преобразование Фурье (Fast Fourier Transform, FFT) находит широкое применение в современной радиотехнике – с его помощью можно не только построить спектр анализируемого сигнала, но оно также необходимо в большинстве операций цифровой обработки сигналов (ЦОС).
Измерение частоты принимаемого сигнала является достаточно сложной задачей для платы Arduino вследствие ее ограниченной вычислительной мощности. Для этой цели частоты используются методы анализа перехода сигнала через ноль – то есть просто измеряется количество переходов сигналов через ноль за определенное время. Но эти методы не работают когда анализируемый сигнал является комбинацией нескольких частот.
В данной статье предпринята попытка написания кода для платы Arduino, позволяющего производить анализ сигналов не вдаваясь глубоко в сущность математических процессов, лежащих в его основе. Этот проект не объясняет принципы быстрого преобразования Фурье (БПФ), но объясняет его использование в Arduino.
Ранее на нашем сайте мы уже рассматривали использование быстрого преобразования Фурье в плате Arduino в следующих проектах:
- анализатор спектра звуковых частот на основе FFT и Arduino;
- 32-полосный анализатор-визуализатор спектра звуковых частот на Arduino.
Но в этих проектах БПФ в Arduino реализовано с помощью специализированных библиотек, в этом же проекте оно реализовано “вручную”, без использования библиотек, и оно позволяет достичь более высоких скоростей БПФ по сравнению с использованием библиотек.
Дискретное преобразование Фурье с высокой скоростью
Чтобы сделать дискретное преобразование Фурье (ДПФ, DFT) более быстрым чем БПФ, James Cooley и John Tukey разработали специальный алгоритм. Некоторые эксперты считают, что этот алгоритм является одним из самых значимых научных прорывов XX века. Данный алгоритм разделяет сигнал на четную и нечетную части, что значительно уменьшает объем необходимых вычислений. При помощи ряда других операций данный алгоритм позволяет свести вычислительную сложность преобразования Фурье к величине NlogN. Обычное ДФП имеет сложность N*N, поэтому данный алгоритм обеспечивает существенное увеличение скорости работы по сравнению с обычным ДПФ при больших N. БПФ также обладает вычислительной сложностью порядка N*logN.
Объяснение данного “революционного” алгоритма ДПФ не является темой настоящей статьи. Посмотреть принципы данного алгоритма вы можете в следующих англоязычных источниках:
https://flylib.com/books/en/2.729.1/derivation_of_the_radix_2_fft_algorithm.html
https://jakevdp.github.io/blog/2013/08/28/understanding-the-fft/
https://cnx.org/contents/8D0YvnW1@7.1:zmcmahhR@7/Decimation-in-time-DIT-Radix-2-FFT
https://en.wikipedia.org/wiki/Fast_Fourier_transform#History
Но наверняка и в русскоязычном сегменте интернета можно найти что-нибудь дельное по этому вопросу. Если кто либо из вас является хорошим специалистом в этой тематике, просьба в комментариях написать сайты, где хорошо написано об этих вопросах.
Общие принципы работы программы
Быстрое преобразование Фурье предполагает многократное вычисление синусов и косинусов. Встроенная в Arduino функция для данных операций не обладает хорошей скоростью работы и требует достаточно много времени для получения приемлемого результата. Чтобы преодолеть недостатки этой функции мы будем значение синуса от 0 до 90 градусов сохранять кратным числу 255, что избавит нас от необходимости использования вещественных чисел (float) и мы сможем хранить эти значения в переменных типа byte, что обеспечивает четырехкратную экономию памяти Arduino. Массив sine_data[ ], который мы будем использовать, необходимо объявить в самом верху программы как глобальный.
Также нам понадобится глобальный массив f_peaks[], который будет обновляться после каждого вызова функции расчета БПФ (FFT). Элемент f_peaks[0] содержит доминирующую (основную) частоту, остальные элементы массива содержат значения в убывающем порядке.
1 2 3 4 5 6 7 8 9 10 11 12 13 |
byte sine_data [91]= { 0, 4, 9, 13, 18, 22, 27, 31, 35, 40, 44, 49, 53, 57, 62, 66, 70, 75, 79, 83, 87, 91, 96, 100, 104, 108, 112, 116, 120, 124, 127, 131, 135, 139, 143, 146, 150, 153, 157, 160, 164, 167, 171, 174, 177, 180, 183, 186, 189, 192, 195, 198, 201, 204, 206, 209, 211, 214, 216, 219, 221, 223, 225, 227, 229, 231, 233, 235, 236, 238, 240, 241, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 253, 254, 254, 254, 255, 255, 255, 255 }; float f_peaks[5]; |
Поскольку мы храним в массиве все значения синуса от 0 до 90 градусов, то мы легко сможем рассчитать любое значение синуса или косинуса. Далее мы будем использовать функцию, которая будет округлять число до первого значения после запятой и возвращать значение из сохраненных данных. Для нашего метода необходима только одна операция с плавающей точкой. От этого недостатка метода тоже можно избавиться если хранить значения синуса, не кратные 255, но это съест много памяти Arduino.
Использованная нами процедура уменьшает точность вычислений, но зато увеличивает их скорость. Для 64 точек она обеспечивает увеличение скорости расчета на 8 мс по сравнению с традиционным подходом, а для 128 точек увеличение скорости расчета составит 20 мс.
Объяснение функции для вычисления БПФ
БПФ может выполняться только на числе отсчетов, кратных 2^n, то есть 4, 8, 16, 32, 64 и т.д. Если число отсчетов не будет равно 2^n, то будут приниматься во внимание только младшие отсчеты. То есть, к примеру, если число отсчетов будет равно 70, то мы будем принимать во внимание только первые 64 отсчета, а остальные отсчеты будут отбрасываться (не учитываться).
Таким образом, БПФ мы будем выполнять на числе отсчетов, равных 2^n, например:
2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048,...
Две переменные вещественного типа out_r и out_im будут занимать достаточно большое количество памяти. Для платы Arduino Nano рассматриваемый алгоритм не будет работать при числе отсчетов большем 128 вследствие недостатка свободной памяти.
1 2 3 4 5 6 7 8 9 10 |
unsigned int data[13]={1,2,4,8,16,32,64,128,256,512,1024,2048}; int a,c1,f,o,x; a=N; for(int i=0;i<12;i++) //calculating the levels { if(data[i]<=a){o=i;} } int in_ps[data[o]]={}; //input for sequencing float out_r[data[o]]={}; //real part of transform float out_im[data[o]]={}; //imaginory part of transform |
Итак, наш алгоритм вычисления БПФ для платы Arduino будет состоять из следующей последовательности действий:
- Код программы генерирует биты в обратном порядке для заданного размера выборки (подробнее о реверсировании битов в шаге 2).
- Входные данные упорядочиваются в соответствии с сгенерированным порядком.
- Выполняется БПФ (FFT).
- Вычисляется амплитуда комплексного числа.
- Обнаруженные максимумы упорядочиваются в убывающем порядке.
- Результаты записываются в массив f_peaks[].
Чтобы получить доступ к другим значениям, находящимся в стороне от пиковой частоты (частоты гармоники), код программы необходимо изменить таким образом, чтобы локальные переменные можно было копировать в некоторые заранее определенные глобальные переменные.
Схема проекта
Для тестирования работы программы можно использовать схему, представленную на следующем рисунке.
Тестирование работы проекта
Для теста написанной программы на ее вход было подано обычно колебание треугольной формы (треугольная волна) с частотой 1.25 Гц. Частота выборки отсчетов для этого колебания составила 10 Гц.
Как показано далее на видео, первоначальные выходные значения программы (raw output) представляют собой рассчитанные значения БПФ, но пока это не те значения, которые нам нужны – они с уменьшенной точностью.
В массиве выходных частот присутствуют частоты 1.25 и 3.75. При этом нет необходимости каждый раз получать точное значение этих частот. Часто эти числа называются элементами разрешения по частоте (frequency bins), то есть выходное значение должно находиться внутри указанных интервалов.
Скорость работы программы для платы Arduino Nano:
16 точек: 4ms
32 точек: 10ms
64 точек: 26ms
128 точек: 53ms
Итак, рассмотренный код программы для вычисления быстрого преобразования Фурье (FFT) можно использовать в приложениях реального времени. Для выполнения расчета необходимо около 30 мс времени. Разрешающая способность программы ограничена числом используемых отсчетов, число которых, в свою очередь, ограничено памятью платы Arduino. Число анализируемых отсчетов и, соответственно, точность работы алгоритма можно улучшить использую платы с большим объемом памяти чем у Arduino Nano - Arduino Mega и т.д.
Исходный код программы (скетча)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 |
/* //Example data: (входные данные, для теста) int data[64]={14, 30, 35, 34, 34, 40, 46, 45, 30, 4, -26, -48, -55, -49, -37, -28, -24, -22, -13, 6, 32, 55, 65, 57, 38, 17, 1, -6, -11, -19, -34, -51, -61, -56, -35, -7, 18, 32, 35, 34, 35, 41, 46, 43, 26, -2, -31, -50, -55, -47, -35, -27, -24, -21, -10, 11, 37, 58, 64, 55, 34, 13, -1, -7}; */ //---------------------------------------------------------------------------// byte sine_data [91]= { 0, 4, 9, 13, 18, 22, 27, 31, 35, 40, 44, 49, 53, 57, 62, 66, 70, 75, 79, 83, 87, 91, 96, 100, 104, 108, 112, 116, 120, 124, 127, 131, 135, 139, 143, 146, 150, 153, 157, 160, 164, 167, 171, 174, 177, 180, 183, 186, 189, 192, 195, //Paste this at top of program (скопируйте это вверх программы 198, 201, 204, 206, 209, 211, 214, 216, 219, 221, 223, 225, 227, 229, 231, 233, 235, 236, 238, 240, 241, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 253, 254, 254, 254, 255, 255, 255, 255 }; float f_peaks[5]; // top 5 frequencies peaks in descending order (первые 5 гармоник в убывающем порядке //---------------------------------------------------------------------------// void setup() { Serial.begin(250000); } void loop() { /* //example FFT(data,64,100); //to get top five value of frequencies of X having 64 sample at 100Hz sampling // для определения первых 5 гармоник при числе отсчетов 64 и частоте отсчетов 100 Гц Serial.println(f_peaks[0]); Serial.println(f_peaks[1]); delay(99999); */ /* after ruing above FFT(), frequencies available at (выходные результаты записываются в) f_peaks[0],f_peaks[1],f_peaks[2],f_peaks[3],f_peaks[4], */ } //-----------------------------FFT Function (функция для вычисления БПФ)----------------------------------------------// float FFT(int in[],int N,float Frequency) { /* Code to perform FFT on arduino, setup: paste sine_data [91] at top of program [global variable], paste FFT function at end of program (скопируйте данные sine_data [91] в верх программы, а функцию FFT в низ программы) Term: 1. in[] : Data array, (массив данных) 2. N : Number of sample (recommended sample size 2,4,8,16,32,64,128...) (число отсчетов) 3. Frequency: sampling frequency required as input (Hz) (частота отсчетов) If sample size is not in power of 2 it will be clipped to lower side of number. i.e, for 150 number of samples, code will consider first 128 sample, remaining sample will be omitted. (если число отсчетов не равно 2 в степени n, по последние отсчеты будут отбрасываться) For Arduino nano, FFT of more than 128 sample not possible due to mamory limitation (64 recomended) (для платы Arduino nano вычисление БПФ для числа отсчетов более 128 невозможно вследствие недостатка памяти, рекомендуется использовать 64 отсчета) For higher Number of sample may arise Mamory related issue, Code by ABHILASH Contact: abhilashpatel121@gmail.com Documentation:https://www.instructables.com/member/abhilash_patel/instructables/ 2/3/2021: change data type of N from float to int for >=256 samples */ unsigned int data[13]={1,2,4,8,16,32,64,128,256,512,1024,2048}; int a,c1,f,o,x; a=N; for(int i=0;i<12;i++) //calculating the levels (рассчитываем уровни) { if(data[i]<=a){o=i;} } int in_ps[data[o]]={}; //input for sequencing (вход последовательности) float out_r[data[o]]={}; //real part of transform (действительная часть преобразования) float out_im[data[o]]={}; //imaginory part of transform (мнимая часть преобразования) x=0; for(int b=0;b<o;b++) // bit reversal (реверсирование битов, смена порядка следования битов) { c1=data[b]; f=data[o]/(c1+c1); for(int j=0;j<c1;j++) { x=x+1; in_ps[x]=in_ps[j]+f; } } for(int i=0;i<data[o];i++) // update input array as per bit reverse order (обновляем входного массив для реверсированного порядка битов) { if(in_ps[i]<a) {out_r[i]=in[in_ps[i]];} if(in_ps[i]>a) {out_r[i]=in[in_ps[i]-a];} } int i10,i11,n1; float e,c,s,tr,ti; for(int i=0;i<o;i++) //fft (БПФ) { i10=data[i]; // overall values of sine/cosine (общее число значений синуса/косинуса) i11=data[o]/data[i+1]; // loop with similar sine cosine(цикл по схожим синусам/косинусам) e=360/data[i+1]; e=0-e; n1=0; for(int j=0;j<i10;j++) { c=cosine(e*j); s=sine(e*j); n1=j; for(int k=0;k<i11;k++) { tr=c*out_r[i10+n1]-s*out_im[i10+n1]; ti=s*out_r[i10+n1]+c*out_im[i10+n1]; out_r[n1+i10]=out_r[n1]-tr; out_r[n1]=out_r[n1]+tr; out_im[n1+i10]=out_im[n1]-ti; out_im[n1]=out_im[n1]+ti; n1=n1+i10+i10; } } } /* for(int i=0;i<data[o];i++) { Serial.print(out_r[i]); Serial.print("\t"); // un comment to print RAW o/p Serial.print(out_im[i]); Serial.println("i"); } */ //---> here onward out_r contains amplitude and our_in conntains frequency (Hz) (out_r содержит амплитуды, а our_in – частоты в Гц) for(int i=0;i<data[o-1];i++) // getting amplitude from compex number (определение амплитуды комплексного числа { out_r[i]=sqrt(out_r[i]*out_r[i]+out_im[i]*out_im[i]); // to increase the speed delete sqrt (для увеличения скорости удалите sqrt) out_im[i]=i*Frequency/N; /* Serial.print(out_im[i]); Serial.print("Hz"); Serial.print("\t"); // un comment to print freuency bin Serial.println(out_r[i]); */ } x=0; // peak detection (обнаружение пиков (гармоник) for(int i=1;i<data[o-1]-1;i++) { if(out_r[i]>out_r[i-1] && out_r[i]>out_r[i+1]) {in_ps[x]=i; //in_ps array used for storage of peak number (массив in_ps используется для хранения номеров гармоник) x=x+1;} } s=0; c=0; for(int i=0;i<x;i++) // re arraange as per magnitude (сортируем в убывающем порядке) { for(int j=c;j<x;j++) { if(out_r[in_ps[i]]<out_r[in_ps[j]]) {s=in_ps[i]; in_ps[i]=in_ps[j]; in_ps[j]=s;} } c=c+1; } for(int i=0;i<5;i++) // updating f_peak array (global variable)with descending order (обновляем массив f_peak (он у нас глобальный) в убывающем порядке { f_peaks[i]=out_im[in_ps[i]]; } } float sine(int i) { int j=i; float out; while(j<0){j=j+360;} while(j>360){j=j-360;} if(j>-1 && j<91){out= sine_data[j];} else if(j>90 && j<181){out= sine_data[180-j];} else if(j>180 && j<271){out= -sine_data[j-180];} else if(j>270 && j<361){out= -sine_data[360-j];} return (out/255); } float cosine(int i) { int j=i; float out; while(j<0){j=j+360;} while(j>360){j=j-360;} if(j>-1 && j<91){out= sine_data[90-j];} else if(j>90 && j<181){out= -sine_data[j-90];} else if(j>180 && j<271){out= -sine_data[270-j];} else if(j>270 && j<361){out= sine_data[j-270];} return (out/255); } //------------------------------------------------------------------------------------// |