FIRフィルタの特性

昨日はMP3音声にFIRフィルタをかけてみました。
その特性をみてみるとこんな感じになります。

で、下のソースのFIRフィルタの部分をこんな感じで奇数番目の符号をかえるようにするとハイパスフィルタになります。

        for(int i = 0; i < data.length - filter.length + 1; ++i){
            for(int j = 0; j < filter.length; ++j){
                proc[i] += data[i + j] * filter[j] * (i % 2 == 0 ? -1 : 1);
            }
        }

特性はこんなん。

ソースはこんな感じで。

import java.awt.*;
import java.awt.geom.GeneralPath;
import java.awt.image.BufferedImage;
import javax.swing.*;
import static java.lang.Math.*;

public class FIR {
    public static void  main(String[] args){
        //処理対象データの作成
        double data[] = new double[100];
        for(int i = 0; i < data.length; ++i){
            data[i] = sin((double)i * i * PI / 2 / data.length);//チャープ信号
        }
        //フィルタ
        double proc[] = new double[data.length];
        //double[] filter = {.2, .2, .2, .2, .2};
        double[] filter = {
            -1.247373789136473e-18,
             1.270350182429102e-02,
            -2.481243022283666e-02,
            -6.381419731491804e-02,
             2.761351394755998e-01,
             6.000000000000000e-01,
             2.761351394755998e-01,
            -6.381419731491804e-02,
            -2.481243022283666e-02,
             1.270350182429102e-02,
            -1.247373789136473e-18}; 
        for(int i = 0; i < data.length - filter.length + 1; ++i){
            for(int j = 0; j < filter.length; ++j){
                proc[i] += data[i + j] * filter[j];
            }
        }
              
        //フーリエ変換
        double[] re = new double[data.length];//実成分
        double[] im = new double[data.length];//虚成分
        double[] pw = new double[data.length];//虚成分
        double max = 0;
        for(int f = 0; f < data.length; ++ f){
            //周期fについてのフーリエ
            re[f] = 0;
            im[f] = 0;
            for(int i = 0; i < data.length; ++i){
                double theta = 2 * PI * f * i / data.length;
                //まちがってた
                //re[f] += proc[i] * cos(theta) + proc[i] * sin(theta);
                //im[f] += proc[i] * cos(theta) - proc[i] * sin(theta);
                //修正 2008/5/1 
                re[f] += proc[i] * cos(theta) - proc[i] * sin(theta);
                im[f] += proc[i] * sin(theta) + proc[i] * cos(theta);//修正 2008/5/17
            }
            re[f] /= data.length;
            im[f] /= data.length;
            pw[f] = sqrt(re[f] * re[f] + im[f] * im[f]);
            if(max < pw[f]) max = pw[f];
        }
        
        //描画
        BufferedImage img = new BufferedImage(400, 300, BufferedImage.TYPE_INT_RGB);
        Graphics2D g = (Graphics2D) img.getGraphics();
        g.setColor(Color.WHITE);
        g.fillRect(0, 0, 400, 300);
        g.setColor(Color.BLACK);
        //源波形
        g.drawString("源波形", 10, 15);
        GeneralPath gp = new GeneralPath();
        gp.moveTo(30, 50);
        for(int i = 0; i < data.length; ++i){
            double x = i * 3 + 30;
            double y = -data[i] * 30 + 50;
            gp.lineTo(x, y);
        }
        g.draw(gp);
        //周波数特性
        g.drawString("周波数特性", 10, 200);
        for(int i = 0; i < re.length / 2; ++i){
            int h = (int) (pw[i] * 85 / max);
            g.fillRect(i * 6 + 30, 290 - h, 4, h + 1);
        }
        //フィルタ波形
        g.drawString("フィルタ波形", 10, 100);
        gp = new GeneralPath();
        gp.moveTo(30, 150);
        for(int i = 0; i < data.length; ++i){
            double x = i * 3 + 30;
            double y = -proc[i] * 30 + 150;
            gp.lineTo(x, y);
        }
        g.draw(gp);
        g.dispose();
        
        //表示
        JFrame f = new JFrame("FIRフィルタ");
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.setSize(400, 350);
        f.add(new JLabel(new ImageIcon(img)));
        f.setVisible(true);
    }
}