フーリエ変換で周波数特性を見る

信号の周波数特性を見るために、フーリエ変換してみます。
一番上が元の波形、その下が周波数特性、で、その下は周波数特性から逆に波形を求めたもので、ここでは元波形に矩形波を入れてます。


ノコギリ波の周波数特性はこんな感じ。

あと、チャープ信号という、広い周波数が一定に入ってる波だとこんな感じです。


フーリエ変換は、関数をsinとcosに分解して微分の計算をやりやすくするためのものなのですが、そのときに周波数ごとの成分が出てくるので、実用的にはもっぱら周波数特性をとってくるのに使われています。
で、コンピュータで処理するので、ここでは離散フーリエ変換(DFT)ということになります。離散フーリエ変換には、効率のよいFFTというアルゴリズムがあるのですが、ここではそれほど速度も求められないし次数も低いし、何より動きがわかりやすいので普通にDFTにしてます。

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

public class Fourier {
    public static void  main(String[] args){
        //処理対象データの作成
        double data[] = new double[100];
        //矩形波
        for(int i = 0; i < data.length; ++i){
            data[i] = (i / 25 % 2 == 0) ? 1 : -1;//矩形波
            //data[i] = i % 50 / 25. - 1;//のこぎり波
            //data[i] = sin((double)i * i * PI / 2 / data.length);//チャープ信号
        }
              
        //フーリエ変換
        double[] re = new double[data.length];//実成分
        double[] im = new double[data.length];//虚成分
        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] += data[i] * cos(theta) + data[i] * sin(theta);
                //im[f] += data[i] * cos(theta) - data[i] * sin(theta);
                //修正 2008/5/1
                re[f] += data[i] * cos(theta) - data[i] * sin(theta);
                im[f] += data[i] * sin(theta) + data[i] * cos(theta);
            }
            re[f] /= data.length;
            im[f] /= data.length;
        }
        
        //振幅スペクトルを得る
        double[] pw = new double[data.length];//振幅スペクトル
        double max = 0;
        for(int f = 0; f < data.length; ++ f){
            pw[f] = sqrt(re[f] * re[f] + im[f] * im[f]);
            if(max < pw[f]) max = pw[f];
        }
        
        //フーリエ逆変換
        double proc[] = new double[data.length];
        for(int f = 0; f < data.length; ++f){
            for(int i = 0; i < data.length; ++i){
                double theta = 2 * PI * f * i / data.length;
                //間違ってた
                //proc[i] += re[f] * cos(theta) - im[f] * sin(theta);
                //修正 2008/5/1
                proc[i] += re[f] * cos(theta) + im[f] * sin(theta);
            }
        }
        
        //描画
        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, 100);
        for(int i = 0; i < re.length / 2; ++i){
            int h = (int) (pw[i] * 85 / max);
            g.fillRect(i * 6 + 30, 190 - h, 4, h + 1);
        }
        //逆変換波形
        g.drawString("フーリエ逆変換", 10, 210);
        gp = new GeneralPath();
        gp.moveTo(30, 250);
        for(int i = 0; i < data.length; ++i){
            double x = i * 3 + 30;
            double y = proc[i] * 30 + 250;
            gp.lineTo(x, y);
        }
        g.draw(gp);
        g.dispose();
        
        //表示
        JFrame f = new JFrame("フーリエ変換");
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.setSize(400, 350);
        f.add(new JLabel(new ImageIcon(img)));
        f.setVisible(true);
    }
}