ヒルベルト変換で瞬間振幅と瞬間周波数

フーリエ変換したやつを、負の周波数の部分を0にして、正の周波数の部分を2倍にするとヒルベルト変換で、これを使うと瞬間振幅とかがとれますよ、と。

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

public class Hilbert {
    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);//チャープ信号
            data[i] = sin(2 * PI * i / data.length * (10 + 10. * i / data.length));
            data[i] *= 1 - sin(PI * i / data.length) * 3 / 4;//振幅の調整
        }

        //フーリエ変換
        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] * sin(theta) + data[i] * cos(theta);
            }
            re[f] /= data.length;
            im[f] /= data.length;
        }
        
        //ヒルベルト変換
        //ω < 0になる成分は0とする(DFTでω < 0の部分はN/2以上の部分に移動している)
        for(int f = data.length / 2; f < data.length; ++f){
            re[f] = 0;
            im[f] = 0;
        }
        //ω > 0の部分を倍にする
        for(int f = 1; f < data.length / 2; ++f){
            re[f] *= 2;
            im[f] *= 2;
        }
        
        //フーリエ逆変換
        double[] proc = new double[data.length];
        double[] hil  = 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);
                hil[i]  += -re[f] * sin(theta) + im[f] * cos(theta);
            }
        }
        
        //瞬時パワー
        double[] power = new double[data.length];
        for(int i = 0; i < data.length; ++i){
            power[i] = proc[i] * proc[i] + hil[i] * hil[i];
        }
        //瞬時周波数
        double[] freq = new double[data.length];
        for(int i = 0; i < data.length - 1; ++i){
            double d = atan2(hil[i + 1], proc[i + 1]) - atan2(hil[i], proc[i]);
            while(d < -PI) d += 2 * PI;
            while(d >= PI) d -= 2 * PI;
            freq[i] = abs(d / (2 * PI));
        }
        
        //描画
        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 / 1; ++i){
            double x = i * 3 + 30;
            double y = -data[i] * 30 + 50;
            gp.lineTo(x, y);
        }
        g.draw(gp);

        //瞬時パワーの描画
        g.drawString("瞬時パワー", 10, 100);
        gp = new GeneralPath();
        gp.moveTo(30, 150);
        for(int i = 0; i < data.length / 1; ++i){
            double x = i * 3 + 30;
            double y = -power[i] * 30 + 180;
            gp.lineTo(x, y);
        }
        g.draw(gp);
        
        //瞬時周波数の描画
        g.drawString("瞬時周波数", 10, 200);
        gp = new GeneralPath();
        gp.moveTo(30, 250);
        for(int i = 0; i < data.length / 1 - 1; ++i){
            double x = i * 3 + 30;
            double y = -freq[i] * 130 + 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);        
    }
}