请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 394967|回复: 2208

[源程序] MATLAB 小波分析源程序(可用) 简要注释

  [复制链接]

新浪微博达人勋

发表于 2012-10-18 22:55:14 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 ccnasyq 于 2013-8-2 16:34 编辑

%最近需要学习小波分析,对论坛的小波分析有关的内容大概看了一下。
%由于缺少注释或者有英文注释(可惜看不懂啊。。),程序运行过程出现问题
%无法找出原因,令人十分困扰。小弟不才,
%借用http://bbs.06climate.com/forum.php?mod=viewthread&tid=7519
%中的程序(此程序经验证可单独使用),简单注释一下,
%以方便后来学习的同志们,不到之处,敬请指出。

clear; %清除变量空间
load q.txt  % 装载q.txt文件到MATLAB
sst = q;    % 获取q.txt文件的句柄,返回给变量sst.
%------------------------------------------------------ Computation  计算
variance = std(sst)^2;  % 计算sst的方差
sst = (sst - mean(sst))/sqrt(variance) ;
n = length(sst);      
dt =1;     %每个Y值之间的时间量 即 x轴的间距,影响图形的拉伸
time = [0:length(sst)-1]*dt + 1963.0 ;  %此处时间“1963”可根据需要自行修改
游客,如果您要查看本帖隐藏内容请回复

% Plot wavelet coefficients of time series
subplot('position',[0.10 0.30 0.80 0.60])
Yticks = 2.^(fix(log2(min(period))):fix(log2(max(period))));
[c,h]=contour(time,log2(period),realpart,'-');  
clabel(c,h);
xlabel('Time (year)')
ylabel('Period (years)')
%title('a) Morlet Wavelet Real Part of the Temperature Anomaly Series')
set(gca,'XLim',xlim(:))
set(gca,'YLim',log2([min(period),max(period)]), ...
    'YDir','default', ...
'YTick',log2(Yticks(:)), ...
'YTickLabel',Yticks)

评分

参与人数 6金钱 +54 贡献 +6 收起 理由
xixinping + 1 很给力!
金玉 + 2 很给力!
italy2114 + 1 很给力!
Aires + 20 + 2
mofangbao + 10 + 2
njzqxt + 20 + 2 谢谢指点。

查看全部评分

本帖被以下淘专辑推荐:

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2014-9-13 12:08:25 | 显示全部楼层
阿吉木柯 2013-6-28 16:34 给LZ的留言:
您好。有个东西想和您讨论一下。
不知道你看过吴洪宝的气候诊断这本书。小波分析这块介绍的很详细。
他也指出了目前很多论文在小波分析的误区。
Morlet小波本来就是复小波,得出来的小波系数是复数,只分析小波系数的实部是不合理的,必须完整的分析使用复小波系数才能恢复信号和表示它的总的变率。所以画出小波系数的实部个人认为是没有意义的,当然墨西哥小波除外,那个是实小波。
而且小波和谐波不一样,如果只分析小波系数模的平方会夸大大尺度成分的强度,所以我看了很多文章也是,在分析了模的平方得出的结果就是,尺度越大的成分强度也很大,这个在气候变化的实际情况其实是不合理的,一般来说,尺度小的振动强度大于尺度大的振动。比如,中纬度地区月平均气温年际变化振幅绝对大于年代际的。
他在文章里也说了(也证明了),小波系数模的平方除以a(伸缩尺度)的平方才能真正反映信号各个周期成分的强度。
我是觉得看了好多好多的小波,碰到的困惑就会越来越多....还是要搞清楚它最终的原理...
就像我知道了要除以a,但是我竟然不知道如何实现,好拙计.....


阿吉木柯提出的问题很值得讨论讨论~
现在有很多人都想学小波分析方法,不考虑其适用性是要不得的。小波分析在大气科学领域的应用已有十几年了,其实这就是一种时频分析工具。比如说一把锤子,适用于钉钉子,你拿去砸核桃就有点不合适了。
阿吉木柯提出的“周期强度”问题,可以采用小波方差检验或小波功率谱检验来解决。用小波分析出的信号变化周期,不经过显著性检验是不可信的。

建议想学习小波分析的童鞋们,先复习一下《统计学》中的周期图、功率谱、交叉谱等知识。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2012-11-15 12:20:44 | 显示全部楼层


>> %最近需要学习小波分析,对论坛的小波分析有关的内容大概看了一下。
%由于缺少注释或者有英文注释(可惜看不懂啊。。),程序运行过程出现问题
%无法找出原因,令人十分困扰。小弟不才,
%借用http://bbs.06climate.com/forum.php?mod=viewthread&tid=7519
%中的程序(此程序经验证可单独使用),简单注释一下,
%以方便后来学习的同志们,不到之处,敬请指出。

clear; %清除变量空间
load q.txt % 装载q.txt文件到MATLAB
sst = q; % 获取q.txt文件的句柄,返回给变量sst.
%------------------------------------------------------ Computation 计算
variance = std(sst)^2; % 计算sst的方差
sst = (sst - mean(sst))/sqrt(variance) ;
n = length(sst);
dt =1; %每个Y值之间的时间量 即 x轴的间距,影响图形的拉伸
time = [0:length(sst)-1]*dt + 1963.0 ; %此处时间“1963”可根据需要自行修改
xlim = [1962,2011]; %此处时间[1962,2011]可根据需要自行修改,即资料时间尺度
pad = 1; % 见matlab小波用法简单翻译http://bbs.06climate.com/forum.p ... &extra=page%3D1
dj = 1/12; %同上
s0 = 2*dt; %同上
j1 =fix((log(n*dt/s0)/log(2))/dj); %同上
mother = 'Morlet'; %同上
% Wavelet transform: 小波变换
[wave,period,scale,coi] = wavelet(sst,dt,pad,dj,s0,j1,mother);%见matlab小波用法简单翻译
power = (abs(wave)).^2 ; % 求小波的功率谱 在此程序里用不上
realpart=real(wave); %求小波的实部
modulus=abs(wave); %求小波的模 在此程序里用不上
phase=atan2(imag(wave),real(wave)); % 求小波的阶 在此程序里用不上
%------------------------------------------------------ Plotting 画图 (以后有机会再详述)
% Plot wavelet coefficients of time series
subplot('position',[0.10 0.30 0.80 0.60])
Yticks = 2.^(fix(log2(min(period))):fix(log2(max(period))));
[c,h]=contour(time,log2(period),realpart,'-');
clabel(c,h);
xlabel('Time (year)')
ylabel('Period (years)')
%title('a) Morlet Wavelet Real Part of the Temperature Anomaly Series')
set(gca,'XLim',xlim(:))
set(gca,'YLim',log2([min(period),max(period)]), ...
'YDir','default', ...
'YTick',log2(Yticks(:)), ...
'YTickLabel',Yticks)
??? Undefined command/function 'wavelet'.

请教下啊,我用这个为什么做了出错啊,能帮忙解答下不
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2012-11-9 14:54:34 | 显示全部楼层
[wave,period,scale,coi] = wavelet(sst,dt,pad,dj,s0,j1,mother);%
Undefined function or method 'wavelet' for input arguments of type 'double'?
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2012-10-19 08:49:21 | 显示全部楼层
本帖最后由 njzqxt 于 2012-10-19 08:50 编辑

图出来了,谢谢!1961-2010年的。
请帮着看看有什么问题。
wa.jpg
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2012-10-23 22:41:41 | 显示全部楼层
请问这个程序用的是每年平均数据还是每月平均数据
密码修改失败请联系微信:mofangbao
回复 支持 1 反对 0

使用道具 举报

新浪微博达人勋

发表于 2012-10-19 08:40:02 | 显示全部楼层
??? mod=viewthread&tid=10545&extra=page%3D1
                      |
Error: The expression to the left of the equals sign is not a valid target for an assignment.
出错啊,加上这句后。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2012-10-19 11:03:28 | 显示全部楼层
支持一下,呵呵
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2012-10-20 08:47:59 | 显示全部楼层
谢谢了,学习一下!!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2012-10-20 08:57:21 | 显示全部楼层
谢谢了,正要学习,这个帮助应该很大!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2012-10-20 13:39:31 | 显示全部楼层
Thank you very mch
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2012-10-28 11:48:07 | 显示全部楼层
好难啊,就是学不会
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2012-10-29 10:00:58 | 显示全部楼层
求翻译下画图部分 实在搞不懂。
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表