チーム15A4
課題名
FFTによる各種フィルタの実装
研究者名
3年 14組 8番 Kazuyuki Kai
3年 14組 9番 Yumeharu Kaji
概要
[PSoC]信号を受け取りフーリエ変換をし、LPFかHPFまたは双方に通し、逆フーリエ変換をし、LCDに再表示する。
[応用版]wav形式ファイルを読み取り、任意のフィルタを通過させたwav形式ファイルを出力する。
PSoCソースコード
#include <stdio.h>
#include <device.h>
#include <math.h>
#define M 32
typedef struct{
double re; double im;
}complex;
void fft(int n, complex a[], complex tmp[]);
void ifft(int n, complex a[], complex tmp[]);
void lpf(int n, complex a[]);
void hpf(int n, complex a[]);
void main()
{
int i; complex a[M], b[M]; LCD_Char_1_Start(); LCD_Char_1_ClearDisplay(); for(i=0;i<M;i++){ b[i].re = 4 * (i % 4); b[i].im = 0; } while(1){
/*LPF*/ if(!LPF_Read()){ fft(M, b, a); lpf(M, b); ifft(M, b, a); while(!LPF_Read()); } /*HPF*/ if(!HPF_Read()){ fft(M, b, a); hpf(M, b); ifft(M, b, a); while(!HPF_Read()); } /*表示*/ for(i=0;i < M;i++){ LCD_Char_1_DrawVerticalBG(1,i,2,b[i].re); } }
}
void fft(int n, complex a[], complex tmp[])
{
int nh,j; double theta; complex x, w;
if(n <= 1){ return; }
theta = -2 * M_PI / n; nh = n / 2;
for(j = 0; j < nh; j++){ tmp[j].re = a[j].re + a[nh+j].re; tmp[j].im = a[j].im + a[nh+j].im; x.re = a[j].re - a[nh+j].re; x.im = a[j].im - a[nh+j].im; w.re = cos(theta*j); w.im = sin(theta*j); tmp[nh+j].re = x.re * w.re - x.im * w.im; tmp[nh+j].im = x.im * w.re + x.re * w.im; }
fft(nh, tmp, a); fft(nh, &tmp[nh], a);
for(j = 0; j < nh; j++){ a[2*j] = tmp[j]; a[2*j+1] = tmp[nh+j]; }
}
void ifft(int n, complex a[], complex tmp[])
{
int j;
for(j = 0; j < n; j++){ a[j].im *= -1; }
fft(n, a, tmp);
for(j = 0; j < n; j++){ a[j].re /= n; a[j].im *= -1; a[j].im /= n; }
}
void lpf(int n, complex a[])
{
int i;
for(i = n * 3 / 4; i < n; i++){ a[i].re = 0; a[i].im = 0; }
}
void hpf(int n, complex a[])
{
int i;
for(i = 0; i < n * 1 / 4; i++){ a[i].re = 0; a[i].im = 0; }
}
応用版ソースコード
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <math.h>
#include <stdint.h>
#define BUF_SIZE 1024
#define N 512 //1フレームのサンプリング数
#define LPF 0
#define HPF 1
#define NF 2
static char buf[BUF_SIZE];
typedef struct{
double re; double im;
}complex;
void fft(int n, complex a[], complex tmp[]);
void ifft(int n, complex a[], complex tmp[]);
void lpf(int n, complex a[]);
void hpf(int n, complex a[]);
void nf(int n, complex a[]);
int main(int argc, char *args[])
{
char buf[BUF_SIZE]; int16_t iBuf; FILE *fp = NULL; FILE *dst = NULL; complex right[N], left[N], tmp[N]; int i, j, cmd;
if(argc != 3){ exit(1); }
if(args[1][0] == '-'){ if(args[1][1] == 'l'){ cmd = LPF; }else if(args[1][1] == 'h'){ cmd = HPF; }else if(args[1][1] == 'n'){ cmd = NF; }else{ exit(1); } }
fp = fopen(args[2], "rb"); assert(fp);
dst = fopen("after.wav", "wb"); assert(dst);
//wavファイルのヘッダ読み取り fread(buf, 52, 1, fp); fwrite(buf, 52, 1, dst);
//フレーム単位で信号処理 while(!feof(fp)){ //元ファイルから1フレーム分読み取り for(i = 0; i < N; i++){ fread(buf, 2, 1, fp); right[i].re = *(int16_t *)buf; right[i].im = 0; fread(buf, 2, 1, fp); left[i].re = *(int16_t *)buf; left[i].im = 0; }
//右チャンネルをFFT->フィルタ->IFFT fft(N, right, tmp); switch(cmd){ case LPF: lpf(N, right); break;
case HPF: hpf(N, right); break;
case NF: nf(N, right); break;
default: break; } ifft(N, right, tmp);
//左チャンネルをFFT->フィルタ->IFFT fft(N, left, tmp); switch(cmd){ case LPF: lpf(N, left); break;
case HPF: hpf(N, left); break;
case NF: nf(N, left); break;
default: break; } ifft(N, left, tmp);
//信号処理後の値を出力するファイルに書き出す for(i = 0; i < N; i++){ iBuf = (int16_t)(right[i].re); fwrite(&iBuf, 2, 1, dst); iBuf = (int16_t)(left[i].re); fwrite(&iBuf, 2, 1, dst); } }
fclose(dst); fclose(fp);
return 0;
}
//フーリエ変換
void fft(int n, complex a[], complex tmp[])
{
int nh,j; double theta; complex x, w;
if(n <= 1){ return; }
theta = -2 * M_PI / n; nh = n / 2;
for(j = 0; j < nh; j++){ tmp[j].re = a[j].re + a[nh+j].re; tmp[j].im = a[j].im + a[nh+j].im; x.re = a[j].re - a[nh+j].re; x.im = a[j].im - a[nh+j].im; w.re = cos(theta*j); w.im = sin(theta*j); tmp[nh+j].re = x.re * w.re - x.im * w.im; tmp[nh+j].im = x.im * w.re + x.re * w.im; }
fft(nh, tmp, a); fft(nh, &tmp[nh], a);
for(j = 0; j < nh; j++){ a[2*j] = tmp[j]; a[2*j+1] = tmp[nh+j]; }
}
//逆フーリエ変換
void ifft(int n, complex a[], complex tmp[])
{
int j;
for(j = 0; j < n; j++){ a[j].im *= -1; }
fft(n, a, tmp);
for(j = 0; j < n; j++){ a[j].re /= n; a[j].im *= -1; a[j].im /= n; }
}
//ローパスフィルタ
void lpf(int n, complex a[])
{
int i;
for(i = n*3/4; i < n; i++){ a[i].re = 0; a[i].im = 0; }
}
//ハイパスフィルタ
void hpf(int n, complex a[])
{
int i;
for(i = 0; i < n/4; i++){ a[i].re = 0; a[i].im = 0; }
}
//ノッチフィルタ(中域周波数のみカット)
void nf(int n, complex a[])
{
int i;
for(i = n/4; i < n*3/4; i++){ a[i].re = 0; a[i].im = 0; }
}
実行方法と仕様
[PSoC]
1.プログラムを実行させる。
2.PSoCのLCDに信号の初期値(初期状態)が表示される。
3.PSoCのボタンを押すとLPFとHPFが実行されるようになっているので、任意でボタンを選択し、押す。(PSoCのLCDを自分側に向けたとき、左がLPF、右がHPF)
4.PSoCのLCDで変化を確認する。
[応用]
1.wav形式ファイルを用意する。(16bitのステレオ音源のみ対応)
2.[コンパイル後のプログラム] [コマンド] [wav形式ファイル]の形式で端末上にて実行させる。
3.コマンドについては"-l"がLPF、"-h"がHPF、"-n"がNFに対応している。また、これ以外のコマンドや上記の形式を満たさない実行については即終了する設定である。
3.フィルタ処理後のwav形式ファイルが生成されるので、視聴または観測等を行う。
各フィルタについては以下のように設定している。
LPF: FFT後の周波数成分のうち上位1/4をすべてカット
HPF: FFT後の周波数成分のうち下位1/4をすべてカット
NF: FFT後の周波数成分のうち上位1/4と下位1/4以外をすべてカット
PSoCの実験結果
図:1.1 初期状態
図:1.2 LPF通過後
図:1.3 HPF通過後
図:1.4 LPF、HPF双方通過させた場合
考察
PSoCのほうのプログラムを実行したところ、LPFとHPFともに成功している(問題がない)と思われるデータが観測された。
また、応用版のほうのプログラムはある音楽データに対して実行した。
それぞれ(LPF、HPF、NF)変化が生じたが一番変化がわかりやすかったのはNFであった。
おそらく、それは今回の実験に使用した音楽データが、視聴する目的で作られたデータであることが関係しているのだと私たちは考えた。
なぜなら、この音楽データはただの波のデータとは違い、人にとって聞きやすい音(すなわち音楽)であり、それは中域の周波数を多く含む(過度な低周波数や高周波数がない)ファイルなのである。
よって、このファイルに対して中盤の周波数を通さないプログラムを実行したところ、逆に聞きにくいものになったと考えた。
(※音楽データについては著作権が発生するためアップロードはしない)
- 最終更新:2015-05-19 10:08:36