チーム12B1
課題名
LPF-FFT解析
研究者名
3年14組4番 大川 拓馬
装置の説明
LPFを通過した信号に対してFFTして、特定の周波数成分を抜き出す装置。
特定周波数の抜出ができればいろいろ応用(例:事象検知器)ができるので、その機能の実装を目指す。
とりあえず基底周波数(D1成分)をLCDに表示する。
装置の構造(ハード部分)
LPFはpsoc creatorのexample projectから流用した。
LPFから出てきた信号を取り込んで波のデータを用意する。
そのためにはサンプリング周波数を用意して、サンプリングする。
精確なサンプリングにはハードウェア割り込みがよいだろうと思い実装を試みたが、
試行錯誤の末に諦めた。無様!
isrが割り込み関係のコンポーネントであることまでは解読できたが、使い方が未解読。時間足りない!
しょうがないのでクロックをピンに出力して、その値でサンプリングを行った(ポーリング)。
ソフトウェアによるものなので他の処理の遅延の影響を受けるため、基底周波数の値が精確である保証はない。
下図の右下に追加してあるのがクロックをピンに出力した配線である。
LPFのフィルタ周波数(cutof周波数)は6kHz
クロックの周波数は128kHz
ポイント数256
だから基底周波数は理想的には
128kHz/256 = 500Hz
となるはずだが、ポーリングなので誤差によりこうはならない。
サンプリング周波数がソフトの処理の遅延の影響を受けて、小さくなるかもしれない。
装置の構造(ソフトウェア部分)
FFT自作した。
再帰呼び出しによるFFTアルゴリズムである。
- 以下fftのソースコード---
#include<stdio.h>
#define _USE_MATH_DEFINES
#include<math.h>
#include<stdlib.h>
#define N 256 //標本化数N
#define exponent(n) log(n)/log(2) //log2(n)を求める
#define depth(n, m) n/m //現在の深さを2のd乗で表現
struct complex_number {
double re; //実数部 double imm;
};
typedef struct complex_number complex; //struct complex_numberではややこしいからcomplexに定義しなおす
complex *w; //回転子を格納する
complex *output; //出力を格納する
int *rev_bit; //ビットリヴァーサルで並べ替えた値を格納する
void initialize(int n) {
int i; for(i = 0;i < n;i++) { output[i].re = 0; output[i].imm = 0; }
}
complex compAdd(complex x, complex y) {
complex result; result.re = x.re + y.re; result.imm = x.imm + y.imm;
return result;
}
complex compMul(complex x, complex y) {
complex result; result.re = x.re * y.re - x.imm * y.imm; result.imm = x.re * y.imm + x.imm * y.re; return result;
}
int *bit_reversal(int n) {
int i,j,x; //配列の動的確保のためにmallocを多用 //外部から標本化数を入力することを想定してプログラムを構成したのでこうなった rev_bit = (int *)malloc(sizeof(int)*n); for(i = 0;i < n;i++) { j = 0; x = i; rev_bit[i] = 0; while(exponent(n) > j) { rev_bit[i] *= 2; if(x&1) { rev_bit[i] += 1; } x /= 2; j++; } } return rev_bit;
}
complex *cal_tf(int n) {
int i; w = (complex *)malloc(sizeof(complex) * n); for(i = 0;i < n;i++) { w[i].re = cos(-2*M_PI*i/n); w[i].imm = sin(-2*M_PI*i/n); } return w;
}
void fft(int n, double input[], int low) {
int i; complex tmp;
//再帰呼び出しの終了条件 if(n == 1) { output[low].re = input[rev_bit[low]]; return; } //再帰呼び出し fft(n/2, input, low); fft(n/2, input, low+n/2); //for(i = 0;i < n;i++) { // printf("output%d:%f+%fi\n", i, output[i].re, output[i].imm); //} //バタフライ演算 for(i = 0;i < n/2;i++) { tmp = output[low+i]; output[low+i] = compAdd(output[low+i], compMul(w[i*depth(N,n)], output[low+i+n/2])); output[low+i+n/2] = compAdd(tmp, compMul(w[(i+n/2)*depth(N,n)], output[low+i+n/2])); }
}
double *amplitude(int n) {
int i; double *amp; amp = (double *)malloc(sizeof(double) * n); for(i = 0;i < n;i++) { amp[i] = sqrt(output[i].re * output[i].re + output[i].imm * output[i].imm); } return amp;
}
double *make_sin(int n) {
int i; double x; double *y; y = (double *)malloc(sizeof(double) * n); for(i = 0;i < n;i++) { x = i/10; y[i] = sin(7 * sin(5 * x)); } return y;
}
void fft_main() {
int i; int *br; double *amp; complex *tf; double *input; //for(i = 0;i < N;i++) { // input[i] = (i < N/2) ? 0 : 1; //} input = make_sin(N); br = bit_reversal(N); for(i = 0;i < N;i++) { printf("%d\n", *(br+i)); } //free(br); tf = cal_tf(N); for(i = 0;i < N;i++) { printf("%f+%fi\n",tf[i].re,tf[i].imm); } //free(tf); output = (complex *)malloc(sizeof(complex) * N); initialize(N); fft(N, input, 0); free(br); free(tf); amp = amplitude(N); for(i = 0;i < N;i++) { printf("output%d:%f+%fi\n", i, output[i].re, output[i].imm); printf("%f\n", amp[i]); } free(output); free(amp); return;
}
int main(void) {
fft_main(); return 0;
}
- ソースコード終わり---
なお、アルゴリズムの動作は確認したものの、その精確さは保証できないので、コピペ流用するときは注意。
だいだいあってる とおもう。
出力outputが1/Nをしていない値なので注意。出力を確認しやすくする(小さくなり過ぎない)ためにこんなことになっている、
使う人は注意だ。
以下にpsocのmain.cのソースコードを記載する。上記のアルゴリズムをどのように組み込んだかがわかると思う。
main関数のループ内がポーリング部分である。
- psocのmain.cのソースコード(一部省略)---
#include <device.h>
#include <VDAC8.h>
#include<math.h>
#include<stdlib.h>
#define N 256
#define exponent(n) log(n)/log(2)
#define depth(n, m) n/m
#define REQUEST_PER_BURST 1u
#define BYTES_PER_BURST 1u
#define UPPER_SRC_ADDRESS CYDEV_PERIPH_BASE
#define UPPER_DEST_ADDRESS CYDEV_PERIPH_BASE
double *fft_main(int n, int input[]);
struct complex_number {
double re; //実数部 double imm;
};
typedef struct complex_number complex; //struct complex_numberではややこしいからcomplexに定義しなおす
complex *w; //回転子を格納する
complex *output; //出力を格納する
int *rev_bit; //ビットリヴァーサルで並べ替えた値を格納する
int input[N];
double *amp_out;
int interrupt_count;
int amp_lcdout;
void DMA_Config(void);
CY_ISR(filterVDAC)
{
/* Convert the 2's complement value to an unsigned 8-bit value * The VDAC expects an unsigned 8-bit value as input. */ VDAC8_Data = (Filter_HOLDAH_REG + 128);
}
void clock_128Hz() {
if(interrupt_count < N) { input[interrupt_count] = Filter_Read24(Filter_CHANNEL_A); interrupt_count++; } else { amp_out = fft_main(N, input); amp_lcdout = amp_out[0]; interrupt_count = 0; LCD_Char_1_Position(1,0); LCD_Char_1_PrintNumber(amp_lcdout); }
}
void initialize(int n) {
int i; for(i = 0;i < n;i++) { output[i].re = 0; output[i].imm = 0; }
}
complex compAdd(complex x, complex y) {
complex result; result.re = x.re + y.re; result.imm = x.imm + y.imm;
return result;
}
complex compMul(complex x, complex y) {
complex result; result.re = x.re * y.re - x.imm * y.imm; result.imm = x.re * y.imm + x.imm * y.re; return result;
}
int *bit_reversal(int n) {
int i,j,x; //配列の動的確保のためにmallocを多用 //外部から標本化数を入力することを想定してプログラムを構成したのでこうなった rev_bit = (int *)malloc(sizeof(int)*n); for(i = 0;i < n;i++) { j = 0; x = i; rev_bit[i] = 0; while(exponent(n) > j) { rev_bit[i] *= 2; if(x&1) { rev_bit[i] += 1; } x /= 2; j++; } } return rev_bit;
}
complex *cal_tf(int n) {
int i; w = (complex *)malloc(sizeof(complex) * n); for(i = 0;i < n;i++) { w[i].re = cos(-2*M_PI*i/n); w[i].imm = sin(-2*M_PI*i/n); } return w;
}
void fft(int n, int input[], int low) {
int i; complex tmp;
//再帰呼び出しの終了条件 if(n == 1) { output[low].re = input[rev_bit[low]]; return; } //再帰呼び出し fft(n/2, input, low); fft(n/2, input, low+n/2); //for(i = 0;i < n;i++) { // printf("output%d:%f+%fi\n", i, output[i].re, output[i].imm); //} //バタフライ演算 for(i = 0;i < n/2;i++) { tmp = output[low+i]; output[low+i] = compAdd(output[low+i], compMul(w[i*depth(N,n)], output[low+i+n/2])); output[low+i+n/2] = compAdd(tmp, compMul(w[(i+n/2)*depth(N,n)], output[low+i+n/2])); }
}
double *amplitude(int n) {
int i; double *amp; amp = (double *)malloc(sizeof(double) * n); for(i = 0;i < n;i++) { amp[i] = sqrt(output[i].re * output[i].re + output[i].imm * output[i].imm); } return amp;
}
double *fft_main(int n, int input[]) {
int i; int *br; double *amp; complex *tf; br = bit_reversal(n); //for(i = 0;i < N;i++) { // printf("%d\n", *(br+i)); //} //free(br); tf = cal_tf(n); //for(i = 0;i < n;i++) { // printf("%f+%fi\n",tf[i].re,tf[i].imm); //} //free(tf); output = (complex *)malloc(sizeof(complex) * n); initialize(n); fft(n, input, 0); free(br); free(tf); amp = amplitude(n); //for(i = 0;i < N;i++) { // printf("output%d:%f+%fi\n", i, output[i].re, output[i].imm); // printf("%f\n", amp[i]); //} //free(output); //free(amp); return amp;
}
void main()
{
/* Start all components used on schematic */ uint8 filter_out; interrupt_count = 0; ADC_DelSig_IRQ_Start(); isr_StartEx(filterVDAC); ADC_DelSig_Start(); ADC_DelSig_StartConvert(); VDAC8_Start(); Opamp_Start(); Filter_Start(); Clock_1_Start(); LCD_Char_1_Start(); LCD_Char_1_ClearDisplay(); LCD_Char_1_Position(0,0); LCD_Char_1_PrintString("D1 spectrum"); LCD_Char_1_Position(1,0); /* User-implemented function to set-up DMA */ DMA_Config(); /* Enable Global Interrupts */ CYGlobalIntEnable;
for(;;) { if(Clock_Out_Read()) { clock_128Hz(); } }
}
- main.cのソースコード終わり---
考察
やはりポーリングではだめだと思う。
基底周波数 = サンプリング周波数/ポイント数
なので、サンプリング周波数で割り込みをかけて標本化が最適だろう。
入力波を変化させてもLCD表示は変わらない。原因は調査中。
だれか割り込みで実装してみてほしい
- 最終更新:2012-12-17 17:41:43