直接上代码,逻辑还算比较清楚
public class SpectrumHelper {
/**
* 计算频谱
* @param data 采样数据
* @param samplingInterval 采样间隔,以秒为单位
* @param needFillZero 是否需要补零
* @param flyLen 参与傅里叶变换的长度
*/
public static SpectrumBean caculateSpectrum(double[] data, double samplingInterval, boolean needFillZero, int flyLen){
if(data == null || data.length == 0){
return null;
}
int dataLen = data.length;
double[] partakeFLYArr = null;//参与傅里叶变换的数组
if(needFillZero){ //需要补零
if(flyLen < dataLen){
return null;
}
partakeFLYArr = new double[flyLen];
System.arraycopy(data, 0, partakeFLYArr, 0, dataLen);
for(int i = dataLen; i < flyLen; i++){
partakeFLYArr[i] = 0;
}
}else{
partakeFLYArr = data;
flyLen = dataLen;
}
//进行傅里叶变换,使用的是commons-math包里面的快速傅里叶算法
FastFourierTransformer fft = new FastFourierTransformer(DftNormalization.STANDARD);
Complex[] complexArr = fft.transform(partakeFLYArr, TransformType.FORWARD);
SpectrumBean spectrumBean = new SpectrumBean();
double[] amplitude = new double[flyLen / 2];
//计算振幅,除以的都是采样点的长度,而不是傅里叶变化的长度
Complex complex = complexArr[0];
amplitude[0] = Math.sqrt(Math.abs(complex.getReal()) * Math.abs(complex.getReal())
+ Math.abs(complex.getImaginary()) * Math.abs(complex.getImaginary())) * 1.0 / dataLen;
for(int i = 1; i < flyLen / 2; i++){
complex = complexArr[i];
amplitude[i] = Math.sqrt(Math.abs(complex.getReal()) * Math.abs(complex.getReal())
+ Math.abs(complex.getImaginary()) * Math.abs(complex.getImaginary())) * 2.0 / dataLen;
}
spectrumBean.setAmplitude(amplitude);//设置振幅
//计算频率分辨率,除以的是傅里叶变换的长度,而不是采样点的长度
double Fs = 1.0 / samplingInterval;//采样频率
spectrumBean.setFrequencyResolution(Fs / flyLen);//设置频率分辨率
return spectrumBean;
}
}
SpectrumBean类表示的是频谱计算的结果,类定义如下:
/**
* 用于保存频谱计算的结果
*/
public class SpectrumBean {
/**
* 频率分辨率
*/
private double frequencyResolution;
/**
* 保存振幅信息,长度是参与傅里叶变换的长度,不是采样点的长度
*/
private double[] amplitude;
public double getFrequencyResolution() {
return frequencyResolution;
}
public void setFrequencyResolution(double frequencyResolution) {
this.frequencyResolution = frequencyResolution;
}
public double[] getAmplitude() {
return amplitude;
}
public void setAmplitude(double[] amplitude) {
this.amplitude = amplitude;
}
}
测试的代码如下:
/**
* 正常的输出波形,正常的傅里叶变换,没有补零的操作
* 频率分辨率 = 采样频率 / 参与傅里叶变换的点数,注意是参与傅里叶变换的点数而非采样点数,因为可能补零
* 傅里叶变换结果的频率计算,第n个点的频率 = 频率分辨率 * (n - 1)
* 所谓频谱,横轴是频率,纵轴是振幅,也就是频率振幅谱
*/
@Test
public void test3(){
//模拟采样
int Fs = 100;//采样频率
int N = 256;//采样点数
double[] wave = generateWaveNormal(Fs, N);
log(Arrays.toString(wave));
//对采样结果做频谱分析
SpectrumBean spectrumBean = SpectrumHelper.caculateSpectrum(wave, 1.0f / Fs, false, 0);
/*for(int i = 0; i < spectrumBean.getAmplitude().length; i++){
log(spectrumBean.getAmplitude()[i]);
}*/
log(spectrumBean.getFrequencyResolution());
}
/**
* 后面补零
* 频率分辨率 = 采样频率 / 参与傅里叶变换的点数,注意是参与傅里叶变换的点数而非采样点数,
* 在补零的情况下,参与傅里叶变换的点数是大于采样点数
* 傅里叶变换结果的频率计算,第n个点的频率 = 频率分辨率 * (n - 1)
* 所谓频谱,横轴是频率,纵轴是振幅,也就是频率振幅谱
*/
@Test
public void test4(){
//模拟采样
int Fs = 100;//采样频率
int N = 256;//采样点数
int WANT_N = 1024;//参与傅里叶变换的点数,实际上只采了256个点,后面的都补零了
double[] wave = generateWaveNormal(Fs, N);
log(Arrays.toString(wave));
//对采样结果做频谱分析
SpectrumBean spectrumBean = SpectrumHelper.caculateSpectrum(wave, 1.0 / Fs, true, WANT_N);
/*for(int i = 0; i < spectrumBean.getAmplitude().length; i++){
log(spectrumBean.getAmplitude()[i]);
}*/
log(spectrumBean.getFrequencyResolution());
}
private void log(Object obj){
System.out.println(obj);
}
/**
* 生成一段模拟信号,模拟采样
* @param fs 采样频率
* @param n 采样点数
* 模拟的这段信号包括的频率有10hz, 15hz, 20hz, 26.5hz,这段信号包括直流信号
* 注意:根据采样定理,采样频率必须大于信号频率的2倍
* 因此,传递过来的采样频率务必要大于26.5 * 2
*/
private double[] generateWaveNormal(float fs, int n) {
double pi = Math.PI;
double[] res = new double[n];
float interval = 1 / fs;//采样间隔,单位是s
for (int i = 0; i < n; i++) {
double t = interval * i;
res[i] = 3 + 1 * Math.cos(2 * pi * 10 * t) + 2 * Math.sin(2 * pi * 15 * t + Math.toRadians(30))
+ 3 * Math.cos(2 * pi * 20 * t + Math.toRadians(-30)) + 4 * Math.sin(2 * pi * 26.5 * t + Math.toRadians(60));
}
return res;
}
测试的结果打印到控制台,然后复制到excel种,使用excel画图
不补零的图如下:
补零的图如下:
结果分析
横坐标频率计算,比如第n个点的频率 = (n -1) * 频率分辨率,频率分辨率的计算见上面
纵坐标表示的是振幅。
频谱图上的信号和原始信号是能正好对应上的。