java版FFT/STFT——转载
2024-05-10 03:05:38
原文资源在此：http://www.ee.columbia.edu/~ronw/code/MEAPsoft/doc/html/files.html
等会我试试效果如何？与我想要的有没有差别。完全公开。良心博文。
00023 package com.meapsoft;
00024
00025
00033 public class FFT {
00034
00035 int n, m;
00036
00037 // Lookup tables. Only need to recompute when size of FFT changes.
00038 double[] cos;
00039 double[] sin;
00040
00041 double[] window;
00042
00043 public FFT(int n) {
00044 this.n = n;
00045 this.m = (int)(Math.log(n) / Math.log(2));
00046
00047 // Make sure n is a power of 2
00048 if(n != (1<<m))
00049 throw new RuntimeException("FFT length must be power of 2");
00050
00051 // precompute tables
00052 cos = new double[n/2];
00053 sin = new double[n/2];
00054
00055 // for(int i=0; i<n/4; i++) {
00056 // cos[i] = Math.cos(-2*Math.PI*i/n);
00057 // sin[n/4-i] = cos[i];
00058 // cos[n/2-i] = -cos[i];
00059 // sin[n/4+i] = cos[i];
00060 // cos[n/2+i] = -cos[i];
00061 // sin[n*3/4-i] = -cos[i];
00062 // cos[n-i] = cos[i];
00063 // sin[n*3/4+i] = -cos[i];
00064 // }
00065
00066 for(int i=0; i<n/2; i++) {
00067 cos[i] = Math.cos(-2*Math.PI*i/n);
00068 sin[i] = Math.sin(-2*Math.PI*i/n);
00069 }
00070
00071 makeWindow();
00072 }
00073
00074 protected void makeWindow() {
00075 // Make a blackman window:
00076 // w(n)=0.42-0.5cos{(2*PI*n)/(N-1)}+0.08cos{(4*PI*n)/(N-1)};
00077 window = new double[n];
00078 for(int i = 0; i < window.length; i++)
00079 window[i] = 0.42 - 0.5 * Math.cos(2*Math.PI*i/(n-1))
00080 + 0.08 * Math.cos(4*Math.PI*i/(n-1));
00081 }
00082
00083 public double[] getWindow() {
00084 return window;
00085 }
00086
00087
00088 /***************************************************************
00089 * fft.c
00090 * Douglas L. Jones
00091 * University of Illinois at Urbana-Champaign
00092 * January 19, 1992
00093 * http://cnx.rice.edu/content/m12016/latest/
00094 *
00095 * fft: in-place radix-2 DIT DFT of a complex input
00096 *
00097 * input:
00098 * n: length of FFT: must be a power of two
00099 * m: n = 2**m
00100 * input/output
00101 * x: double array of length n with real part of data
00102 * y: double array of length n with imag part of data
00103 *
00104 * Permission to copy and use this program is granted
00105 * as long as this header is included.
00106 ****************************************************************/
00107 public void fft(double[] x, double[] y)
00108 {
00109 int i,j,k,n1,n2,a;
00110 double c,s,e,t1,t2;
00111
00112
00113 // Bit-reverse
00114 j = 0;
00115 n2 = n/2;
00116 for (i=1; i < n - 1; i++) {
00117 n1 = n2;
00118 while ( j >= n1 ) {
00119 j = j - n1;
00120 n1 = n1/2;
00121 }
00122 j = j + n1;
00123
00124 if (i < j) {
00125 t1 = x[i];
00126 x[i] = x[j];
00127 x[j] = t1;
00128 t1 = y[i];
00129 y[i] = y[j];
00130 y[j] = t1;
00131 }
00132 }
00133
00134 // FFT
00135 n1 = 0;
00136 n2 = 1;
00137
00138 for (i=0; i < m; i++) {
00139 n1 = n2;
00140 n2 = n2 + n2;
00141 a = 0;
00142
00143 for (j=0; j < n1; j++) {
00144 c = cos[a];
00145 s = sin[a];
00146 a += 1 << (m-i-1);
00147
00148 for (k=j; k < n; k=k+n2) {
00149 t1 = c*x[k+n1] - s*y[k+n1];
00150 t2 = s*x[k+n1] + c*y[k+n1];
00151 x[k+n1] = x[k] - t1;
00152 y[k+n1] = y[k] - t2;
00153 x[k] = x[k] + t1;
00154 y[k] = y[k] + t2;
00155 }
00156 }
00157 }
00158 }
00159
00160
00161
00162
00163 // Test the FFT to make sure it's working
00164 public static void main(String[] args) {
00165 int N = 8;
00166
00167 FFT fft = new FFT(N);
00168
00169 double[] window = fft.getWindow();
00170 double[] re = new double[N];
00171 double[] im = new double[N];
00172
00173 // Impulse
00174 re[0] = 1; im[0] = 0;
00175 for(int i=1; i<N; i++)
00176 re[i] = im[i] = 0;
00177 beforeAfter(fft, re, im);
00178
00179 // Nyquist
00180 for(int i=0; i<N; i++) {
00181 re[i] = Math.pow(-1, i);
00182 im[i] = 0;
00183 }
00184 beforeAfter(fft, re, im);
00185
00186 // Single sin
00187 for(int i=0; i<N; i++) {
00188 re[i] = Math.cos(2*Math.PI*i / N);
00189 im[i] = 0;
00190 }
00191 beforeAfter(fft, re, im);
00192
00193 // Ramp
00194 for(int i=0; i<N; i++) {
00195 re[i] = i;
00196 im[i] = 0;
00197 }
00198 beforeAfter(fft, re, im);
00199
00200 long time = System.currentTimeMillis();
00201 double iter = 30000;
00202 for(int i=0; i<iter; i++)
00203 fft.fft(re,im);
00204 time = System.currentTimeMillis() - time;
00205 System.out.println("Averaged " + (time/iter) + "ms per iteration");
00206 }
00207
00208 protected static void beforeAfter(FFT fft, double[] re, double[] im) {
00209 System.out.println("Before: ");
00210 printReIm(re, im);
00211 fft.fft(re, im);
00212 System.out.println("After: ");
00213 printReIm(re, im);
00214 }
00215
00216 protected static void printReIm(double[] re, double[] im) {
00217 System.out.print("Re: [");
00218 for(int i=0; i<re.length; i++)
00219 System.out.print(((int)(re[i]*1000)/1000.0) + " ");
00220
00221 System.out.print("]\nIm: [");
00222 for(int i=0; i<im.length; i++)
00223 System.out.print(((int)(im[i]*1000)/1000.0) + " ");
00224
00225 System.out.println("]");
00226 }
00227 }STFT在此：
00023 package com.meapsoft;
00024
00025 import java.io.IOException;
00026 import java.util.ArrayList;
00027 import java.util.Arrays;
00028
00029 import javax.sound.sampled.AudioFormat;
00030 import javax.sound.sampled.AudioInputStream;
00031
00039 public class STFT {
00040 AudioInputStream input;
00041 AudioFormat format;
00042 int bytesPerWavFrame, frameLen;
00043 ArrayList listeners = new ArrayList();
00044 double[] re, im, window;
00045 static double log10 = Math.log(10);
00046 static double epsilon = 1e-9; // avoid log of zero
00047 RingMatrix freq, time;
00048 FFT fft;
00049
00050 static double rmsTarget = 0.08;
00051 static double rmsAlpha = 0.001;
00052 double rms = 1;
00053
00054 public float samplingRate;
00055 public int nhop;
00056
00057 // The line should be open, but not started yet.
00058 public STFT(AudioInputStream input, int frameLen, int hopsize, int history) {
00059 this(input, frameLen, history);
00060
00061 nhop = hopsize;
00062 }
00063
00064 // The line should be open, but not started yet.
00065 public STFT(AudioInputStream input, int frameLen, int history) {
00066 freq = new RingMatrix(frameLen/2+1, history);
00067 time = new RingMatrix(frameLen, history);
00068 this.frameLen = frameLen;
00069
00070 this.input = input;
00071 format = input.getFormat();
00072 bytesPerWavFrame = format.getFrameSize();
00073
00074 samplingRate = format.getSampleRate();
00075
00076 fft = new FFT(frameLen);
00077
00078 this.re = new double[frameLen];
00079 this.im = new double[frameLen];
00080 this.window = fft.getWindow();
00081 for(int i=0; i<im.length; i++)
00082 im[i] = 0;
00083
00084 nhop = frameLen;
00085
00086 samplingRate = input.getFormat().getSampleRate();
00087 }
00088
00089 // returns the total number of bytes read
00090 public long start() {
00091 byte[] b = new byte[bytesPerWavFrame * frameLen];
00092 Arrays.fill(b, (byte)0);
00093
00094 int noverlapBytes = (frameLen-nhop)*bytesPerWavFrame;
00095 int nhopBytes = nhop*bytesPerWavFrame;
00096
00097 int totalBytesRead = 0;
00098 int bytesRead = 22;
00099 while(bytesRead > 0) {
00100 if(nhop > 0)
00101 {
00102 // shift b so our overlap works out
00103 for(int x = 0; x < noverlapBytes; x++)
00104 b[x] = b[x+nhopBytes];
00105 }
00106
00107 try {
00108 bytesRead = input.read(b, noverlapBytes, nhopBytes);
00109 totalBytesRead += bytesRead;
00110 } catch(IOException ioe) {
00111 ioe.printStackTrace();
00112 return totalBytesRead;
00113 }
00114
00115 // store the unwindowed waveform for getSamples function
00116 double[] wav = time.checkOutColumn();
00117 bytes2doubles(b, wav);
00118
00119 // Normalize rms using a moving average estimate of it
00120 // Calculate current rms
00121 double rmsCur = 0;
00122 for(int i=0; i<wav.length; i++)
00123 rmsCur += wav[i]*wav[i];
00124 rmsCur = Math.sqrt(rmsCur / wav.length);
00125
00126 // update moving average
00127 rms = rmsAlpha*rmsCur + (1-rmsAlpha)*rms;
00128
00129 // normalize by rms
00130 for(int i=0; i<wav.length; i++)
00131 wav[i] = wav[i] * rmsTarget / rms;
00132
00133 time.checkInColumn(wav);
00134
00135 // window waveform
00136 for(int i=0; i<wav.length; i++)
00137 re[i] = window[i] * wav[i];
00138
00139 // take fft
00140 fft.fft(re, im);
00141
00142 // Calculate magnitude
00143 double[] mag = freq.checkOutColumn();
00144 for(int i=0; i<mag.length; i++)
00145 mag[i] = 10*Math.log(re[i]*re[i] + im[i]*im[i] + epsilon) / log10;
00146
00147 // clear im[]
00148 Arrays.fill(im, 0);
00149
00150 // Tell everyone concerned that we've added another frame
00151 long frAddr = freq.checkInColumn(mag);
00152 notifyListeners(frAddr);
00153 }
00154
00155 // let the frame listeners know that we're done reading:
00156 notifyListeners(-1);
00157
00158 return totalBytesRead;
00159 }
00160
00161
00162 // Get the waveform samples from frames frStart to frEnd-1
00163 public double[] getSamples(long frStart, long frEnd) {
00164 long sampStart = fr2Samp(frStart);
00165 long sampEnd = fr2Samp(frEnd);
00166
00167 double[] x = new double[(int)(sampEnd - sampStart)];
00168
00169 for(int fr=0; fr < frEnd-frStart; fr++) {
00170 double[] frame = time.getColumn(frStart+fr);
00171 if(frame == null) continue;
00172 // only the first nhop samples of frame are valid
00173 for(int i=0; i<nhop; i++)
00174 x[(int)(fr2Samp(fr+frStart)-sampStart + i)] = frame[i];
00175 }
00176
00177 return x;
00178 }
00179
00180 // Convert an address in frames into an address in samples
00181 public long fr2Samp(long frAddr) {
00182 return nhop * frAddr;
00183 }
00184
00185 // Convert an address in samples into an address in frames
00186 public long samp2fr(long sampAddr) {
00187 return sampAddr/nhop;
00188 }
00189
00190 public double[] getFrame(long frAddr) { return freq.getColumn(frAddr); }
00191 public void setFrame(long frAddr, double[] dat) { freq.setColumn(frAddr, dat); }
00192 public int getColumns() { return freq.getColumns(); }
00193 public int getRows() { return freq.getRows(); }
00194
00195 // Dealing with FrameListeners
00196 public void addFrameListener(FrameListener fl) {
00197 listeners.add(fl);
00198 }
00199 public void removeFrameListener(FrameListener fl) {
00200 listeners.remove(fl);
00201 }
00202 public void notifyListeners(long frAddr) {
00203 for(int i=0; i<listeners.size(); i++) {
00204 FrameListener list = (FrameListener) listeners.get(i);
00205 list.newFrame(this, frAddr);
00206 }
00207 }
00208
00209
00210 // Convert a byte stream into a stream of doubles. If it's stereo,
00211 // the channels will be interleaved with each other in the double
00212 // stream, as in the byte stream.
00213 public void bytes2doubles(byte[] audioBytes, double[] audioData) {
00214 if (format.getSampleSizeInBits() == 16) {
00215 if (format.isBigEndian()) {
00216 for (int i = 0; i < audioData.length; i++) {
00217 /* First byte is MSB (high order) */
00218 int MSB = (int) audioBytes[2*i];
00219 /* Second byte is LSB (low order) */
00220 int LSB = (int) audioBytes[2*i+1];
00221 audioData[i] = ((double)(MSB << 8 | (255 & LSB)))
00222 / 32768.0;
00223 }
00224 } else {
00225 for (int i = 0; i < audioData.length; i++) {
00226 /* First byte is LSB (low order) */
00227 int LSB = (int) audioBytes[2*i];
00228 /* Second byte is MSB (high order) */
00229 int MSB = (int) audioBytes[2*i+1];
00230 audioData[i] = ((double)(MSB << 8 | (255 & LSB)))
00231 / 32768.0;
00232 }
00233 }
00234 } else if (format.getSampleSizeInBits() == 8) {
00235 int nlengthInSamples = audioBytes.length;
00236 if (format.getEncoding().toString().startsWith("PCM_SIGN")) {
00237 for (int i = 0; i < audioBytes.length; i++) {
00238 audioData[i] = audioBytes[i] / 128.0;
00239 }
00240 } else {
00241 for (int i = 0; i < audioBytes.length; i++) {
00242 audioData[i] = (audioBytes[i] - 128) / 128.0;
00243 }
00244 }
00245 }
00246 }
00247
00248 // Convert an address in frames to an address in seconds
00249 public double fr2Seconds(long frAddr)
00250 {
00251 return(fr2Samp(frAddr)/samplingRate);
00252 }
00253
00254 // Convert an address in seconds to an address in frames
00255 public long seconds2fr(double sec)
00256 {
00257 return(samp2fr((long)(sec*samplingRate)));
00258 }
00259
00260
00264 public static RingMatrix getSTFT(double[] samples, int nfft)
00265 {
00266 return STFT.getSTFT(samples, nfft, nfft);
00267 }
00268
00272 public static RingMatrix getSTFT(double[] samples, int nfft, int nhop)
00273 {
00274 RingMatrix freq = new RingMatrix(nfft/2+1, samples.length/nhop);
00275
00276 FFT fft = new FFT(nfft);
00277 double[] window = fft.getWindow();
00278
00279 double[] wav = new double[nfft];
00280 double rms = 1;
00281 for(int currFrame = 0; currFrame < samples.length/nhop; currFrame++)
00282 {
00283 // zero pad if we run out of samples:
00284 int zeroPadLen = currFrame*nhop + wav.length - samples.length;
00285 if(zeroPadLen < 0)
00286 zeroPadLen = 0;
00287 int wavLen = wav.length - zeroPadLen;
00288
00289 //for(int i = 0; i<wav.length; i++)
00290 // wav[i] = samples[currFrame*nhop + i];
00291 for(int i = 0; i < wavLen; i++)
00292 wav[i] = samples[currFrame*nhop + i];
00293 for(int i = wavLen; i < wav.length; i++)
00294 wav[i] = 0;
00295
00296 // Normalize rms using a moving average estimate of it
00297 // Calculate current rms
00298 double rmsCur = 0;
00299 for(int i=0; i<wav.length; i++)
00300 rmsCur += wav[i]*wav[i];
00301 rmsCur = Math.sqrt(rmsCur / wav.length);
00302
00303 // update moving average
00304 rms = rmsAlpha*rmsCur + (1-rmsAlpha)*rms;
00305
00306 // normalize by rms
00307 for(int i=0; i<wav.length; i++)
00308 wav[i] = wav[i] * rmsTarget / rms;
00309
00310 // window waveform
00311 double[] re = new double[wav.length];
00312 double[] im = new double[wav.length];
00313 for(int i=0; i<wav.length; i++)
00314 {
00315 re[i] = window[i] * wav[i];
00316 im[i] = 0;
00317 }
00318
00319 // take fft
00320 fft.fft(re, im);
00321
00322 // Calculate magnitude
00323 double[] mag = freq.checkOutColumn();
00324 for(int i=0; i<mag.length; i++)
00325 mag[i] = 10*Math.log(re[i]*re[i] + im[i]*im[i] + epsilon) / log10;
00326
00327 freq.checkInColumn(mag);
00328 }
00329
00330 return freq;
00331 }
00332
00333
00342 public int readFrames(long nframes)
00343 {
00344 byte[] b = new byte[bytesPerWavFrame * frameLen];
00345 Arrays.fill(b, (byte)0);
00346
00347 int noverlapBytes = (frameLen-nhop)*bytesPerWavFrame;
00348 int nhopBytes = nhop*bytesPerWavFrame;
00349
00350 int bytesRead = 22;
00351 int nFramesRead = 0;
00352 while(nFramesRead <= nframes)
00353 {
00354 if(nhop > 0)
00355 {
00356 // shift b so our overlap works out
00357 for(int x = 0; x < noverlapBytes; x++)
00358 b[x] = b[x+nhopBytes];
00359 }
00360
00361 try
00362 {
00363 input.read(b, noverlapBytes, nhopBytes);
00364 nFramesRead++;
00365 }
00366 catch(IOException ioe)
00367 {
00368 ioe.printStackTrace();
00369 return nFramesRead;
00370 }
00371
00372 // store the unwindowed waveform for getSamples function
00373 double[] wav = time.checkOutColumn();
00374 bytes2doubles(b, wav);
00375
00376 // Normalize rms using a moving average estimate of it
00377 // Calculate current rms
00378 double rmsCur = 0;
00379 for(int i=0; i<wav.length; i++)
00380 rmsCur += wav[i]*wav[i];
00381 rmsCur = Math.sqrt(rmsCur / wav.length);
00382
00383 // update moving average
00384 rms = rmsAlpha*rmsCur + (1-rmsAlpha)*rms;
00385
00386 // normalize by rms
00387 for(int i=0; i<wav.length; i++)
00388 wav[i] = wav[i] * rmsTarget / rms;
00389
00390 time.checkInColumn(wav);
00391
00392 // window waveform
00393 for(int i=0; i<wav.length; i++)
00394 re[i] = window[i] * wav[i];
00395
00396 // take fft
00397 fft.fft(re, im);
00398
00399 // Calculate magnitude
00400 double[] mag = freq.checkOutColumn();
00401 for(int i=0; i<mag.length; i++)
00402 mag[i] = 10*Math.log(re[i]*re[i] + im[i]*im[i] + epsilon) / log10;
00403
00404 // clear im[]
00405 Arrays.fill(im, 0);
00406
00407 // Tell everyone concerned that we've added another frame
00408 long frAddr = freq.checkInColumn(mag);
00409 notifyListeners(frAddr);
00410 }
00411
00412 return nFramesRead;
00413 }
00414
00418 public long getLastFrameAddress()
00419 {
00420 return freq.nextFrAddr-1;
00421 }
00422
00423 public void stop() throws IOException
00424 {
00425 input.close();
00426 }
00427 }另外有相关问题可以加入QQ群讨论，不设微信群
QQ群：868373192
语音深度学习群
java版FFT/STFT——转载相关推荐
- 服务端工程师入门与进阶 Java 版
前言 欢迎加入我们.这是一份针对实习生/毕业生的服务端开发入门与进阶指南.遇到问题及时问你的 mentor 或者直接问我. 建议: 尽量用google查找技术资料. 有问题在stackoverflow ...
- base64转码java版
base64转码java版 package com.net.util;import java.io.FileInputStream; import java.io.FileOutputStream; ...
- Java版 QQ空间自动登录无需拷贝cookie一天抓取30WQQ说说数据流程分析【转】
Java版 QQ空间自动登录无需拷贝cookie一天抓取30WQQ说说数据&流程分析 QQ空间说说抓取难度比较大,花了一个星期才研究清楚! 代码请移步到GitHub GitHub地址:http ...
- 线索二叉树原理及前序、中序线索化（Java版）
转载 原文地址:https://blog.csdn.net/UncleMing5371/article/details/54176252 一.线索二叉树原理 前面介绍二叉树原理及特殊二叉树文章中提到, ...
- 【java项目实践】具体解释Ajax工作原理以及实现异步验证username是否存在+源代码下载(java版)...
一年前,从不知道Ajax是什么,伴随着不断的积累,到如今常常使用,逐渐有了深入的认识. 今天,假设想开发一个更加人性化,友好,无刷新,交互性更强的网页,那您的目标一定是Ajax. 介绍 在具体讨论Aj ...
- JAVA版StarDict星际译王简单实现
由胡正开发的星际译王是Linux平台上很强大的一个开源的翻译软件(也有Windows版本的)支持多种词库.多种语言版本.尤其词库设计比较合理.之前看到一篇博文<星际译王词库应用-自制英汉词典&g ...
- Java版世界时钟示例
Java版世界时钟示例 这是一个Java版的世界时钟示例,移植自Gerrit创建的同名Swing应用(http://www.jug-muenster.de/swing-worldclock-427 ) ...
- 经典十大排序算法（含升序降序，基数排序含负数排序）【Java版完整代码】【建议收藏系列】
经典十大排序算法[Java版完整代码] 写在前面的话 十大排序算法对比 冒泡排序 快速排序 直接选择排序 堆排序 归并排序 插入排序 希尔排序 计数排序 桶排序 基数排序 完整测试类 写在前面的话 ...
- 20165234 [第二届构建之法论坛] 预培训文档(Java版) 学习总结
[第二届构建之法论坛] 预培训文档(Java版) 学习总结 我通读并学习了此文档,并且动手实践了一遍.以下是我学习过程的记录~ Part1.配置环境 配置JDK 原文中提到了2个容易被混淆的概念 JD ...
最新文章
- 万万没想到 I 这 7 件超酷的事情，让开发更有效率
- 箱线图怎么判断异常值_原创【六西格玛工具解读】02——箱线图(Boxplot)
- redis的观察者模式----------发布订阅功能
- 八大排序算法的 Python 实现
- Redis介绍及部署在CentOS7上（一） 1
- 所有特征在不同分类之间、 train和test之间的列分布差异(图形绘制)
- 第十节：进一步扩展两种安全校验方式
- linux mysql提示1045_linux mysql ERROR 1045
- 数据结构之线索二叉树
- DSO 中的Windowed Optimization
- idea之springboot端口被占用/跳转到login
- keras的训练引擎：train_array.py和train_generator.py
- CCF NOI1041 志愿者选拔
- MySQL 中文的乱码问题
- android下最强的3款pdf阅读器测评
- 01_测试基础知识---功能测试常用方法/正交表的使用
- 【重识云原生】第三章云存储第一节——分布式云存储总述
- 音视频中的语音信号处理技术
- 怎样把自己喜欢的微信表情包(动态)导出来，我三岁半的表弟都会...
- 国内硕士申请加拿大计算机博士难度,加拿大硕士和博士真的那么难申请吗？