フーリエ変換したやつを、負の周波数の部分を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); } }