注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

阿弥陀佛

街树飘影未见尘 潭月潜水了无声 般若观照心空静...

 
 
 

日志

 
 
关于我

一直从事气象预报、服务建模实践应用。 注重气象物理场、实况场、地理信息、本体知识库、分布式气象内容管理系统建立。 对Barnes客观分析, 小波,计算神经网络、信任传播、贝叶斯推理、专家系统、网络本体语言有一定体会。 一直使用Java、Delphi、Prolog、SQL编程。

网易考拉推荐

CWT.java  

2013-02-25 18:36:17|  分类: 小波 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
/* From: "A SURVEY OF COMPUTATIONAL PHYSICS" 
   by RH Landau, MJ Paez, and CC BORDEIANU 
   Copyright Princeton University Press, Princeton, 2008.
   Electronic Materials copyright: R Landau, Oregon State Univ, 2008;
   MJ Paez, Univ Antioquia, 2008; and CC BORDEIANU, Univ Bucharest, 2008.
   Support by National Science Foundation                               
   */
// CWT.java  Continuous Wavelet Transform. 
import java.io.*;
import java.util.*;

public class CWT  {
    
  public static void main(String[] argv) throws IOException, FileNotFoundException {
    int i, j, k, o, t, n=120, m=100; 
    double c[][]= new double[m][n];  double input[] = new double[1024], omega1 = 1.;
    double tau, dtau, omega, domega, x, WTreal, WTimag, max, omega2= 5.,tau1= -81.92;
    double PsiReal[]= new double[16834], PsiImag[]= new double[16834]; 
  
    dtau = 1./m;
    System.out.println("Transform in CWT.dat, signal in CWTin.dat");  
    PrintWriter w = new PrintWriter (new FileOutputStream("CWT.dat"), true);
    PrintWriter q = new PrintWriter (new FileOutputStream("CWTin.dat"), true);
    for( t=0; t < 750; t++ ) {
      if( t > 0    && t <= 250 ) input[t] =  5*Math.sin(6.28*t);
      if( t >= 250 && t <= 750 ) input[t] = 10*Math.sin(2.*6.28*t);
      q.println(""+input[t]+"");
    } 
    // Psi(t) = Psi((t-tau)/s) = Psi((t-tau)*omega), tau2 = tau1 + n*dtau  
    domega = Math.pow(omega2/omega1, 1.0/m);                               // Scaling
    omega = omega1;
    for ( i = 0; i < m; i++ )  {                          // Compute daughter wavelet
     tau = tau1;
      for (o=0; o<16834; o++)  {                // For signals up to 2^13 = 8192 long
        PsiReal[o] = WaveletReal( tau*omega );
        PsiImag[o] = WaveletImag( tau*omega );
        tau = tau + dtau;                                              // Translation
      }            
      for ( j=0; j<n; j++ ) {                                          // Compute CWT
        WTreal = 0.;
        WTimag = 0.;
        for (o=0;o<input.length;o++)  { 
          WTreal += input[o]*PsiReal[8192-(j*input.length)/n+o];
          WTimag += input[o]*PsiImag[8192-(j*input.length)/n+o]; 
        }
        c[i][j] = Math.sqrt(WTreal*WTreal+WTimag*WTimag);
      }
      omega = omega*domega;                                                // Scaling
    }
    max = 0.0001;       
    for (i = 0; i<m; i++) {
      for ( j=0; j<n; j++ ) {                                          // Renormalize
        if ( c[i][j] > max )   max = c[i][j];                
        w.println(" "+c[i][j]/max+" "); 
      }      
      w.println(""); 
    } 
  }
    
  public static double WaveletReal( double t) {                  // Re Morlet wavelet
    double sigma = 4.0;
    return Math.cos(2.0*Math.PI*t)* Math.exp(-t*t/(2.*sigma*sigma))
                                        / (sigma*Math.sqrt(2.*Math.PI));
  }
        
  public static double WaveletImag( double t ) {                 // Im Morlet wavelet
    double sigma = 4.0;
    return Math.sin(2.0*Math.PI*t)* Math.exp(-1.*t*t/(2.*sigma*sigma))
                                       / (sigma*Math.sqrt(2.0*Math.PI));
  }     
}
  评论这张
 
阅读(397)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017