hhh-wavelet feature5
- 格式:pdf
- 大小:149.43 KB
- 文档页数:4
[收稿日期]2021-12-19 [修回日期]2023-04-19[基金项目]蚌埠医学院自然科学研究重点项目(2021byzd073)[作者单位]蚌埠医学院第一附属医院1.超声科,3.放射科,4.肿瘤外科,安徽蚌埠233004;2.蚌埠医学院超声医学教研室,安徽蚌埠233030[作者简介]苏 蕾(1987-),女,硕士,主治医师.[通信作者]孙医学,主任医师,副教授.E⁃mail:133****1158@163.com[文章编号]1000⁃2200(2023)08⁃1101⁃04㊃影像医学㊃超声影像组学对乳腺BI⁃RADS 3类及以上结节良恶性鉴别诊断的应用价值苏 蕾1,2,周牧野3,李 阳1,郭 婕1,侯迎迎1,马小开4,孙医学1[摘要]目的:探讨超声影像组学对乳腺BI⁃RADS 3类及以上结节良恶性的鉴别诊断㊂方法:回顾性分析经超声诊断为BI⁃RADS 3类及以上结节164个,并获得其病理结果,将病灶划分为训练集(n =123)㊁测试集(n =41)㊂基于超声影像组学中的特征提取工程提取图像中感兴趣区的高维特征参数并利用特征降维筛选可靠特征建立模型进行分类研究,训练集用于建立模型,测试集用于验证模型准确性㊂结果:超声组学所建立的模型鉴别乳腺结节良恶性的敏感性㊁特异性㊁准确性及ROC 曲线下面积分别为87.51%㊁84.62%㊁85.00%㊁0.854(0.707~0.946)㊂结论:超声影像组学在鉴别乳腺BI⁃RADS 3类及以上结节的良恶性上具有良好的诊断价值㊂[关键词]超声影像组学;乳腺结节;BI⁃RADS;鉴别诊断[中图法分类号]R 445.1;R 737.9 [文献标志码]A DOI :10.13898/ki.issn.1000⁃2200.2023.08.019Application value of ultrasound radiomics in differential diagnosis of benign and malignant breast BI⁃RADS category 3and above nodulesSU Lei 1,2,ZHOU Mu⁃ye 3,LI Yang 1,GUO Jie 1,HOU Ying⁃ying 1,MA Xiao⁃kai 4,SUN Yi⁃xue 1(1.Department of Ultrasound ,3.Department of Radiology ,4.Department of Urgical Oncology ,The First Affiliated Hospital of Bengbu Medical College ,Bengbu Anhui 233004;2.Ultrasound Medicine Teachingand Research Office ,Bengbu Medical College ,Bengbu Anhui 233030,China )[Abstract ]Objective :To evaluate the differential diagnosis of benign and malignant breast BI⁃RADS category 3and above nodules by ultrasound radiomics.Methods :A total of 164breast nodules of BI⁃RADS category 3and above diagnosed by ultrasound were retrospectively analyzed.The lesions were divided into training set (n =123)and test set (n =41).Base on feature extraction engineering of ultrasound radiomics,the high latitude feature parameters of the region of interest in the image were extracted,and the feature dimension reduction was used to screen the reliable features which established the model for classification research.The training set was used to establish the model,and the test set was used to verify the accuracy of the model.Results :The sensitivity,specificity,accuracy and area under ROC curve of the model established by ultrasound radiomics were 87.51%,84.62%,85.00%and 0.854(0.707-0.946),respectively.Conclusions :Ultrasound radiomics has certain diagnostic value in differentiating benign and malignant breast BI⁃RADS category 3and above nodules.[Key words ]ultrasound radiomics;breast nodule;BI⁃RADS;differential diagnosis 近年来,乳腺癌的发病年龄呈年轻化趋势,对其治疗的方案也更倾向个体化㊁精准化,因此,在保证病人生存获益的前提下,同时要满足女性的正常生理㊁生活需求,乳腺癌的手术范围不断缩小,手术方式更加微创,以期最大程度保留病人乳房的功能,减少术后可能带来的并发症,这就对乳腺癌的精确诊断提出更大的挑战㊂目前,对于乳腺疾病的筛查普遍利用超声来进行,并依据BI⁃RADS 分类[1]标准进行诊断,但不同医生的诊断存在一定的主观差异性㊂而影像组学作为新型影像图像技术,通过提取㊁量化图像的高维特征,对定量数据进行分析,避免主观误判㊂本研究基于超声影像组学,通过构建统计学模型来分析及预测超声诊断为BI⁃RADS 3类及以上乳腺结节的良恶性㊂现作报道㊂1 资料与方法1.1 一般资料 选取2018年6月至2021年11月蚌埠医学院第一附属医院收治的乳腺结节病人164例,均为女性,年龄26~68岁,平均(41.7±8.4)岁,共164个病灶,且经超声诊断均为BI⁃RADS 3类及3类以上,后经病理证实良性病灶71例,恶性病灶93例㊂纳入标准:(1)所有病人图像清晰,肿瘤形态能完全显示在图像中;(2)病人行超声检查获取图像前均未行任何辅助治疗;(3)病人除乳腺肿瘤外均未有其他肿瘤㊂排除标准:(1)获取病理前提示存在远处转移或合并其他肿瘤;(2)超声图像存在大量声影,肿瘤较大导致感兴趣区(ROI)不能完全勾画㊂本项目经蚌埠医学院伦理委员会批准,伦科批字[2022]第186号㊂1.2 超声检查及方法 选用GE S8彩色超声诊断仪,高频线阵探头11L⁃D,频率3~11MHz,对病人进行图像采集㊂病人选取仰卧位,行双侧乳腺常规超声扫查,记录图像中结节的大小㊁数目及位置,根据病灶形态㊁大小㊁边界㊁内部回声㊁后方回声衰减特征㊁钙化及病灶与周围组织之间关系等进行BI⁃RADS分类,多枚结节共存的情况下,只保留分类等级最高的结节作为ROI选取,所有图像经由2名具有7年以上乳腺超声诊断经验的医生采集并评估,最终综合一致性意见得到分类结果㊂采集图像时调节频率㊁焦点㊁增益及时间补偿曲线使图像达到最佳成像质量㊂每个乳腺肿瘤病人选取最大长轴切面图像,导出为DICOM格式文件㊂1.3 图像预处理 原始图像由960mm×720mm×1mm体素尺寸(Voxel Size)㊁1mm×1mm×1mm体素间距(Voxel Spcaing)和[0,255]的灰度值构成㊂在特征提取之前对每幅图像都进行相同的标准化预处理,将原始3D图像转换为960mm×720mm像素尺寸,1mm×1mm像素间距,且对灰度值进行Z⁃score标准化处理,获得均值为0㊁标准差为1的灰度图像㊂转换公式:(x⁃μ)/σ,x为原始像素点值,μ为原始像素均值,σ为原始像素标准差㊂1.4 ROI勾画㊁特征提取及数据预处理 所有图像由一名副高职称医生使用3D⁃slicer(Version 4.11.0,)软件对病灶ROI进行手动勾画,并由另一名副高职称医生核查,导出肿瘤图像及ROI掩模(见图1)㊂使用Python(Version3.7.6, )的Pyradiomics程序库[2]对所有原始图像进行小波变换(Wavelet filter)及指数变换(Exponential filter),获得16幅变换图像,同时提取原始图像㊁小波变换及指数变换图像的2D组学特征㊂对所有提取的超声组学特征值进行Z⁃score标准化处理㊂1.5 特征建模与统计分析 1.5.1 特征初步筛选 为了消除特征之间多重共线性对预测结果的影响,组内采取Pearson相关系数对重复性较高的特征进行筛选,相关系数r>0.9的特征认为与其他特征高度线性相关,并予以剔除㊂组间采用t检验和Mmann⁃whitney检验,筛选P>0.05的特征予以剔除,剩余特征入组待下一步降维㊁筛选㊂1.5.2 特征降维及建模 将164个病灶按照3∶1的比例依据BI⁃RADS分类随机分层分为训练集(n=123)和测试集(n=41)㊂训练集用于特征降维㊁筛选及建立模型,测试集用于评估模型㊂运用L1正则化(LASSO)算法对剩余入组特征进行压缩㊁降维,采用十折交叉验证筛选最优lambda值,特征系数不为0的参数作为最终纳入超声影像组学评分模型Radscore的变量㊂1.6 统计学方法 采用Mmann⁃whitney检验㊁t检验㊁χ2检验和Pearson相关性分析,运用错配率㊁ROC曲线及曲线下面积(AUC)评价二元logistic回归模型的预测能力㊂2 结果2.1 一般资料 164个病灶中,共71例良性病灶, 93例恶性病灶㊂训练集123例,测试集41例㊂训练集与测试集组结节最大径㊁良恶性分布以及类别分布差异均无统计学意义(P>0.05)(见表1)㊂表1 训练集和测试集良恶性病灶数及BI⁃RADS分类结果分组n结节最大径[M(P25,P75)]/mm结节良恶性良性 恶性 BI⁃RADS分类 3类 4类 5类 训练集12325(6,35)5370307419测试集 4123(5,36)182310247合计16425(6,36)7193409826χ2 0.78*0.010.07P >0.05>0.05>0.05 *示Mmann⁃whitney检验2.2 特征筛选及模型构建 在训练组上使用特征提取器将每幅图像提取567个特征,组内运用Pearson相关系数剔除r>0.9的特征,共有289个特征予以保留,组间经过t检验和Mmann⁃whitney检验,筛选出P<0.05的特征㊂组内和组间初筛后剩余97个特征㊂采用LASSO进一步对以上特征进行降维,将特征数压缩成8个影像组学特征和一个截距值,同时获得最优Lambda值0.03328,错配率0.0724(见图2),8个特征包括1个一阶特征㊁2个2D形体学特征及5个纹理特征,5个原始图像特征,1个对数变换特征,2个小波变换特征,见表2中①~⑧㊂ 利用上述特征建立二元logistic回归方程, Radscore(Radscore:1.Radscore=1.29203907×①-1.82709327×②+0.17149944×③+0.016 01477×④+0.47144247×⑤-0.17895320×⑥-0.14634798×⑦+0.34389966×⑧-⑨)为最终保留的特征与其对应系数乘积之和再与截距(Intercept)相加,最终预测变量为是否为乳腺癌概率P(二分类逻辑回归模型:概率P=1÷(1+ e-Radscore),特征系数为数据标准化过后的回归系数及截距值(见表2)㊂表2 特征筛选表序号特征名称特征系数①original_shape2D_Elongation 1.29203907②original_shape2D_Sphericity-1.82709327③original_firstorder_Entropy0.17149944④original_glcm_ClusterProminence0.01601477⑤original_glrlm_RunEntropy0.47144247⑥exponential_glrlm_RunLengthNonUniformity-0.17895320⑦Wavelet-LH_gldm_DependenceVariance-0.14634798⑧Wavelet-LH_glszm_ZoneEntropy0.34389966⑨Intercept-0.968663222.3 模型的诊断效能 在测试集41例病人超声图像上对影像组学模型进行验证㊂模型的敏感性㊁特异性㊁阴性预测值㊁阳性预测值㊁准确性㊁Youden指数分别为87.51%㊁84.62%㊁90.91%㊁77.78%㊁85.00%㊁0.7033,AUC为0.854(0.707~0.946) (P<0.01)(见图3)㊂结果表明超声影像组学鉴别良恶性乳腺病变具有良好的诊断效能㊂3 讨论 目前,对乳腺疾病的影像学诊断包括超声㊁钼靶和MRI等㊂MRI虽然具有多序列成像特点,但对于较小病灶的检测具有一定的局限性,而且需要依赖Gd增强和弥散加权成像来辅助判断㊂钼靶X线具有高分辨率成像特点,但为重叠的静态图像,当乳腺腺体增多时,X线发现病灶难度增加,无法进行精准诊断,且具有辐射性㊂而超声具有图像清晰无重叠,检查方便快速,无辐射,发现病灶准确等特点,常作为首选的检查手段㊂超声影像组学作为新的影像组学技术,近些年被逐渐重视㊂超声基于声波成像,不同于其他影像学的密度及氢质子成像方式,运用组学对超声图像分析方面可以获得其他影像组学不同的特征信息,从而构建不同的模型,因此超声影像组学是乳腺疾病诊断技术中不可或缺的一部分㊂本研究采用超声影像组学对超声图像进行特征提取和特征降维㊁分析,提取特征分为训练集和测试集数据㊂在训练集中动态的降低模型分类结果与真实分类结果之间的误差从而筛选Lamdba值和最佳特征参数,构建逻辑回归模型,最终因变量返回为乳腺癌的概率P㊂测试集上运用ROC曲线及AUC方法验证逻辑回归模型分类的准确性㊂模型分类的敏感性和特异性较高,分别为87.51%和84.62%㊂筛选出的8个超声组学特征用来构建模型,提取了包括一阶直方图特征㊁2D形态学特征和高阶纹理特征[3]㊂其中original_shape2D_Elongation和original_ shape2D_Sphericity为形态学特征,Elongation将图像分成两个主体部分,计算其伸长率,Sphericity为反映图像球形率的特征㊂original_firstorder_Entropy为一阶直方图特征熵,反映图像灰阶分布的随机性㊂而其余5个均是高阶纹理特征[4],original_glcm_ ClusterProminence反映ROI集群突出的特征,纹理集群越突兀,不对称性越大,其值越大㊂original_ glrlm_RunEntropy反映灰度行程的不确定性,其值增加,ROI的异质性和不均匀性增加㊂exponential_ glrlm_RunLengthNonUniformity测量灰度游程相似性,较低值表示灰度长程更加均匀㊂Wavelet⁃LH_ gldm_DependenceVariance测量ROI区域灰度依赖方差大小㊂Wavelet⁃LH_glszm_ZoneEntropy反映ROI灰度分布区域的不确定性和随机性,值越大表明随机性越大㊂相对于一阶直方图特征与形态学特征,纹理特征作为诊断㊁鉴别良恶性肿瘤具有明显的意义,肿瘤纹理特征表示肿瘤图像的像素排布与空间像素关系[5],代表肿瘤生长发展的纹理异质性表现[6]㊂以上这些表现在图像上无法辨别,同时也难以发现这些微观信息与疾病发生的内在联系㊂超声组学对鉴别乳腺病变的良恶性的研究较少,多数集中在对甲状腺结节的分析[7]㊂本研究数据纳入了BI⁃RADS3类及以上的结节,因3类及以上结节均具有恶性潜能,而BI⁃RADS分级高也并非完全是恶性,有研究[8]表明,BI⁃RADS4类结节良性率高达62%,恶性可能性为2%~95%[9]㊂大部分影像组学对鉴别乳腺癌良恶性的研究集中在MRI 及X线,既往资料显示MRI及X线利用影像组学鉴别乳腺结节的AUC分别为0.75~0.81[10]㊁0.76~ 0.86[11],而本研究的AUC为0.85,更是高于MRI 组学,与X线组学基本接近㊂影像组学利用机器学习思想建立模型预测乳腺疾病的性质已经得到证实[12-13],本研究通过对乳腺BI⁃RADS3类及以上结节的分析发现,超声影像组学作为影像组学对乳腺肿瘤良恶性鉴别诊断具有重要价值,同时拥有快速㊁经济以及无辐射等特点,但研究存在一定的局限性㊂首先,本次研究选取样本为单中心样本,需要多中心样本来增加模型的泛化性;其次,影像图像不同于自然图像容易获得,同时由于病灶的异质性等原因不能运用自然图像的图像增强技术来处理,存在一定的模型过拟合问题,今后将开展多中心㊁大样本研究,同时改进超声影像图像的提取技术,减少过拟合等问题,增加模型预测的泛化性及准确性㊂[参考文献][1] EGHTEDARI M,CHONG A,RAKOW⁃PENNER R,et al.Currentstatus and future of BI⁃RADS in multimodality imaging,from theAJR special series on radiology reporting and data systems[J].AJR Am J Roentgenol,2021,216(4):860.[2] GRIETHUYSEN J,FEDOROV A,PARMAR C,et al.Computational radiomics system to decode the radiographicphenotype[J].Cancer Res,2017,77(21):e104. [3] CHEN Q,XIA J,ZHANG J.Identify the triple⁃negative and non⁃triple⁃negative breast cancer by using texture features of medicaleultrasonic image:a STROBE⁃compliant study[J].Medicine(Baltimore),2021,100(22):e25878.[4] ALOBAIDLI S,MCQUAID S,SOUTH C,et al.The role of textureanalysis in imaging as an outcome predictor and potential tool inradiotherapy treatment planning[J].Br J Radiol,2014,87(1042):20140369.[5] ZHANG B,SONG L,YIN J.Texture analysis of DCE⁃MRIintratumoral subregions to identify benign and malignant breasttumors[J].Front Oncol,2021,11:688182.[6] DAMASCELLI A,GALLIVANONE F,CRISTEL G,et al.Advanced imaging analysis in prostate MRI:building a radiomicsignature to predict tumor aggressiveness[J].Diagnostics(Basel),2021,11(4):594.[7] 胡小玲,冉海涛.超声影像组学评估甲状腺乳头状癌颈部淋巴结转移[J].中国超声医学杂志,2022,38(4):4. [8] ELVERICI E,BARCA AN,AKTAS H,et al.Nonpalpable BI⁃RADS4breast lesions:sonographic findings and pathologycorrelation[J].Diagn Interv Radiol,2015,21(3):189. [9] 吴佳辉,鹿皎,卓晓英,等.声触诊弹性成像技术联合超声造影优化乳腺BI⁃RADS4类病灶诊断的价值研究[J].徐州医科大学学报,2023,43(3):196.[10] GIBBS P,ONISHI N,SADINSKI M,et al.Characterization of sub⁃1cm breast lesions using radiomics analysis[J].J Magn ResonImaging,2019,50(5):1468.[11] LI Z,YU L,WANG X,et al.Diagnostic performance ofmammographic texture analysis in the differential diagnosis ofbenign and malignant breast tumors[J].Clin Breast Cancer,2018,18(4):e621.[12] 王瑛,陈英格,叶素敏,等.超声影像组学标签预测三阴性乳腺癌的价值[J].中国医学影像学杂志,2021,29(1):35. [13] 徐敏,马宜传,张书海,等.基于双模态MRI影像组学术前预测浸润性乳腺癌腋窝淋巴结转移[J].蚌埠医学院学报,2021,46(12):1763.(本文编辑 周洋)。
\documentclass{article}\usepackage{graphicx}\usepackage[round]{natbib}\bibliographystyle{plainnat}\usepackage[pdfstartview=FitH,%bookmarksnumbered=true,bookmarksopen=true,%colorlinks=true,pdfborder=001,citecolor=blue,%linkcolor=blue,urlcolor=blue]{hyperref}\begin{document}\title{Research plan under the Post-doctorate program at xx University}%\subtitle{aa}\author{Robert He}\date{2008/04/23}\maketitle\section{Research Title}~~~~Crustal seismic anisotropy in the xx using Moho P-to-S converted phases.\section{Research Background \& Purposes}~~~~Shear-wave splitting analyses provide us a new way to study the seismic structure and mantle dynamics in the crust and mantle. The crustal anisotropy is developed due to various reasons including lattice-preferred orientation (LPO) of mineral crystals and oriented cracks.\newlineTraditionally, the earthquakes occurring in the curst and the subducting plates are selected to determine the seismic anisotropy of the crust. However, none of these methods can help us to assess the anisotropy in the whole crust. Because crustal earthquakes mostly are located in the upper crust, they do not provide information of lower crust. On the other hand, earthquakes in the subducting plates provide information of the whole crust but combined with upper mantle. However, it’s difficult to extract the sole contri bution of the crust from the measurement. Fortunately P-to-S converted waves (Ps) at the Moho are ideal for investigation of crustal seismic anisotropy since they are influenced only by the medium above the Moho.Moho. Figure \ref{crustalspliting}~schematically shows the effects of shear wave splitting on Moho Ps phases. Initially, a near-vertically incident P wave generates a radially polarized converted shear wave at the crust-mantle boundary. The phases, polarized into fast and slow directions, progressively split in time as they propagate through the anisotropic media. Here, the Ps waves can be obtained from teleseismic receiver function analysis.%%\begin{figure}[htbp]\begin{center}\includegraphics[width=0.47\textwidth]{crustalsplit.png}\caption{The effects of shear wave splitting in the Moho P to S converted phase. Top shows a schematic seismogram in the fast/slow coordinate system with split horizontal Ps components.(cited from: McNamara and Owens, 1993)}\label{crustalspliting}\end{center}\end{figure}%%The Korean Peninsula is composed of three major Precambrian massifs, the Nangrim, Gyeongii, and Yeongnam massifs(Fig.\ref{geomap}). The Pyeongbuk-Gaema Massif forms the southern part of Liao-Gaema Massif of southern Manchuria, and the Gyeonggi and Mt. Sobaeksan massifs of the peninsula are correlated with the Shandong and Fujian Massifs of China.%\begin{figure}[htbp]\begin{center}\includegraphics[width=0.755\textwidth]{geo.png}\caption{Simplified geologic map. NCB: North China block; SCB: South China block.(cited from: Choi et al., 2006)}\label{geomap}\end{center}\end{figure}%Our purpose of the study is to measure the shear wave splitting parameters in the crust of the Korean Peninsula. The shear wave splitting parameters include the splitting time of shear energybetween the fast and slow directions, as well as fast-axis azimuthal direction in the Korean Peninsula. These two parameters provide us constraints on the mechanism causing the crustal anisotropy. From the splitting time, the layer thickness of anisotropy will be estimated. Whether crustal anisotropy mainly contributed by upper or lower crustal or both will be determined. Based on the fast-axis azimuthal direction, the tectonic relation between northeastern China and the Korean peninsula will be discussed.\section{Research Methods}~~~~Several methods have been introduced for calculation of receiver functions. An iterative deconvolution technique may be useful for this study since it produces more stable receiver function results than others. The foundation of the iterative deconvolution approach is aleast-squares minimization of the difference between the observed horizontal seismogram and a predicted signal generated by the convolution of an iteratively updated spike train with the vertical-component seismogram. First, the vertical component is cross-correlated with the radial component to estimate the lag of the first and largest spike in the receiver function (the optimal time is that of the largest peak in the absolute sense in the cross-correlation signal). Then the convolution of the current estimate of the receiver function with the vertical-component seismogram is subtracted from the radial-component seismogram, and the procedure is repeated to estimate other spike lags and amplitudes. With each additional spike in the receiver function, the misfit between the vertical and receiver-function convolution and the radial component seismogram is reduced, and the iteration halts when the reduction in misfit with additional spikes becomes insignificant.\newlineFor all measurement methods of shear-wave splitting, time window of waveform should be selected. Conventionally the shear-wave analysis window is picked manually. However, manual window selection is laborious and also very subjective; in many cases different windows give very different results.\newlineIn our study, the automated S-wave splitting technique will be used, which improves the quality of shear-wave splitting measurement and removes the subjectivity in window selection. First, the splitting analysis is performed for a range of window lengths. Then a cluster analysis isapplied in order to find the window range in which the measurements are stable. Once clusters of stable results are found, the measurement with the lowest error in the cluster with the lowest variance is presented for the analysis result.\section{Expected results \& their contributions}~~~~First, the teleseismic receiver functions(RFs) of all stations including radial and transverse RFs can be gained. Based on the analysis of RFs, the crustal thickness can be estimated in the Korean Peninsula. Then most of the expected results are the shear-wave splitting parameters from RFs analysis in the crust beneath the Korean Peninsula. The thickness of anisotropic layer will be estimated in the region when the observed anisotropy is assumed from a layer of lower crustal material.All the results will help us to understand the crustal anisotropy source.\newlineCrustal anisotropy can be interpreted as an indicator of the crustal stress/strain regime. In addition, since SKS splitting can offer the anisotropy information contributed by the upper mantle but combined with the crust, the sole anisotropy of the upper mantle can be attracted from the measurement of SKS splitting based on the crustal splitting result.%\cite{frogge2007}%%%\citep{frogge2008}%%%\citep{s-frogge2007}% 5. References\begin{thebibliography}{99}\item Burdick, L. J. and C. A. Langston, 1977, Modeling crustal structure through the use of converted phases in teleseismic body waveforms, \textit{Bull. Seismol. Soc. Am.}, 67:677-691.\item Cho, H-M. et al., 2006, Crustal velocity structure across the southern Korean Peninsula from seismic refraction survey, \textit{Geophy. Res. Lett.} 33, doi:10.1029/2005GL025145.\item Cho, K. H. et al., 2007, Imaging the upper crust of the Korean peninsula by surface-wave tomography, \textit{Bull. Seismol. Soc. Am.}, 97:198-207.\item Choi, S. et al., 2006, Tectonic relation between northeastern China and the Korean peninsula revealed by interpretation of GRACE satellite gravity data, \textit{Gondwana Research}, 9:62-67.\item Chough, S. K. et al., 2000, Tectonic and sedimentary evolution of the Korean peninsula: a review and new view, \textit{Earth-Science Reviews}, 52:175-235.\item Crampin, S., 1981, A review of wave motion in anisotropic and cracked elastic-medium, \textit{Wave Motion}, 3:343-391.\item Fouch, M. J. and S. Rondenay, 2006, Seismic anisotropy beneath stable continental interiors, \textit{Phys. Earth Planet. Int.}, 158:292-320.\item Herquel, G. et al., 1995, Anisotropy and crustal thickness of Northern-Tibet. New constraints for tectonic modeling, \textit{Geophys. Res. Lett.}, 22(14):1 925-1 928.\item Iidaka, T. and F. Niu, 2001, Mantle and crust anisotropy in the eastern China region inferred from waveform splitting of SKS and PpSms, \textit{Earth Planets Space}, 53:159-168.\item Kaneshima, S., 1990, Original of crustal anisotropy: Shear wave splitting studies in Japan, \textit{J. Geophys. Res.}, 95:11 121-11 133.\item Kim, K. et al., 2007, Crustal structure of the Southern Korean Peninsula from seismic wave generated by large explosions in 2002 and 2004, \textit{Pure appl. Geophys.}, 164:97-113.\item Kosarev, G. L. et al., 1984, Anisotropy of the mantle inferred from observations of P to S converted waves, \textit{Geophys. J. Roy. Astron. Soc.}, 76:209-220.\item Levin, V. and J. Park, 1997, Crustal anisotropy in the Ural Mountains foredeep from teleseismic receiver functions, \textit{Geophys. Res. Lett.}, 24(11):1 283 1286.\item Ligorria, J. P. and C. J. Ammon, 1995, Iterative deconvolution and receiver-function estimation. \textit{Bull. Seismol. Soc. Am.}, 89:1 395-1 400.\item Mcnamara, D. E. and T. J. Owens, 1993, Azimuthal shear wave velocity anisotropy in the basin and range province using Moho Ps converted phases, \textit{J. Geophys. Res.}, 98:12003-12 017.\item Peng, X. and E. D. Humphreys, 1997, Moho dip and crustal anisotropy in northwestern Nevada from teleseismic receiver functions, \textit{Bull. Seismol. Soc. Am.}, 87(3):745-754.\item Sadidkhouy, A. et al., 2006, Crustal seismic anisotropy in the south-central Alborz region using Moho Ps converted phases, \textit{J. Earth \& Space Physics}, 32(3):23-32.\item Silver, P. G. and W. W. Chan, 1991, Shear wave splitting and subcontinental mantle deformation, \textit{J. Geophys. Res.},96:16 429-16454.\item Teanby, N. A. et al., 2004, Automation of shear wave splitting measurement using cluster analysis, \textit{Bull. Seismol. Soc. Am.}, 94:453-463.\item Vinnik, L. and J-P. Montagner, 1996, Shear wave splitting in the mantle Ps phases,\textit{Geophys. Res. Lett.}, 23(18):2 449- 2 452.\item Yoo, H. J. et al., 2007, Imaging the three-dimensional crust of the Korean peninsula by joint inversion of surface-wave dispersion of teleseismic receiver functions, \textit{Bull. Seismol. Soc. Am.}, 97(3):1 002-1 011.\item Zhu, L., and H. Kanamori, 2000, Moho depth variation in Southern California from teleseismic receiver functions, \textit{J. Geophys. Res.}, doi :10.1029/1999JB900322, 105:2 969-2 980.%%%%\end{document}。
Users ManUal Model #93711EnglishCongratulations on your purchase of the Celestron NexImage 5 Solar System imaging camera. Your neximage camera comes with the following:+ NexImage 5 Camera + 1.25” nose piece +CD including NexImage iCap and RegiStax processing software + USB CableRecommended Minimum system requirements:+ Pentium IV , 2.0 GHz, 1GB RAM + Graphics card with 24 or 32 bit +Windows XP , Windows Vista, Windows 7 (32 & 64 bit)+ DirectX 9.0c or higher1.25” Nose Piece USB Cable NexImage 5 Camera NexImage iCap and Registax processing softwareFig 11Quick start1. Place the CD into your computer’s CD rom drive.2. Install both the NexImage iCap and RegiStax softwareonto your computer.3. Thread the 1.25” nose piece into the body of the NexImage camera. see Fig 2.4. Slide the 1.25” nose piece of the camera into the eyepiece barrel of your telescope. see Fig 3.5. Plug the small end of the USB cable into the port on the back of the camera.6. Plug camera into the USB port of your computer. Once connected, the Found New Hardwaremessage will appear. Follow the Installation Wizard until the device has been successfully installed.2Fig 2Fig 3Capturing images1. Double-click the NexImage iCap icon on your computer’s desktop to start the program. see Fig 4.2. If the NexImage 5 camera is not already detected, select the camera and press OK. see Fig 5.3. Select an easy target like the moon to begin with. Center and focus your telescope on the specific feature you wish to image.4.You should be able to see light displayed in the Preview window of the iCap software.3Fig 4Fig 55. Use the telescopes focuser to focus the image until object is visible and sharp. see Fig6.6. Use the Gain (A ) and Exposure (B ) settings to make sure that the image is not under or over exposed.7. Select the video format (C ). For the highest resolution, select 2592x1944.8. Select the frames per second (FPS) for your video. The FPS box (D ) will always display the maximum frames per second for the selected resolution setting.9. Press the Video File button (E ) to bring up the Recording Setting box.see Fig 7 and 8.4Fig 6Fig 7Fig 8B D C10. Press the Video File (E ) to select the destination location for your video.11. Click on the Advanced tab. You can either select the amount of time or number of frames you wish to record. see Fig 9.12. Press the Record button (F ) when ready to begin recording video. see Fig 10.13. Once recording is complete, press the check video button (g ) to view thecapture image.See the Help file (H) to learn more about the advanced features of iCap software.Now that your video is captured, you are ready to process it into one high resolution image.5Fig 9Fig 10Processing images1. Double-click the RegiStax icon on your computer’s desktop to start the program. see Fig 11.2. Press the Select button (A ) and select the video that was just recorded. see Fig 12.3. Press the Set Alignment Points button (B ) to have RegiStax automatically select the alignment positions over your image.4. Press the ALIGN button (C ) to begin the alignment process.5. Select the Best Frame option under the Limit Setup box (D ) and enter a number that represents one-half of the number of image frames that were captured, i.e. if 300 frames were captured then enter 150. Press the LIMIT button (E).6Fig 11Fig 126. You will automatically advance to the STACK screen. Accept the default settings and press the STACK button (F ). see Fig 13.77. Next moved to the wavelet processing screen by clicking on the Wavelet tab (g ). see Fig 14.8. On the wavelet page use the wavelet-sliders (h ) to enhancestored in the image whereas the larger numbered wavelets to reveal the desired amount of detail for your image.9. Fig 13FGTo explore the many other features of Registax and view helpful tutorials, go to the Registax homepage at: http://www.astronomie.be/registax/index.htmlhow neximage worksNexImage utilizes a light sensitive imaging sensor to capture streaming video of any solar system object. This video can easily be viewed as hundreds of individual images (frames) that can be digitally stacked to significantly reduce the electric “noise” inherent in video chips and bring out the unseen fine detail (signal) hidden within your image. The NexImage camera takes advantage of the fact that the signal to noise ratio of your stacked composite image is proportional to the square root of the number of frames combined. This means that stacking as few as 16 frames will reduce the grainy noise of the composite image by 4 times. While stacking as many 900 frames will improve the image by 30 times! However, stacking the individual frames is only half the power of the NexImage imager. With the included software package, each individual frame is analyzed for quality to filter out those frames most affected (blurred) by poor atmospheric “seeing”. This form of after-the-fact adaptive optics, leaves only the sharpest, clearest frames to be stacked and aligned into a high quality image. Finally, powerful processing features automatically break the image up into individual unsharp mask layers that can be used to bring out tremendous detail and reveal final images that will rival those taken with astronomical CCD cameras costing thousands of dollars.The BasicsFocusingAs with all astrophotography, sharp focus is essential for high quality results. Although there are many techniques and devices for focusing your telescope, the human eye still remains one of the best detectors of subtle changes in detail. One advantage that video imaging has over imaging with more sophisticated (and expensive) CCD cameras is the speed in which it can display its image. Focusing NexImage is more similar to focusing an eyepiece than a CCD camera. Unlike with long exposure cameras you don’t have to wait many seconds to see the effect of a focus change.8Focusing TipsTo achieve best focus, concentrate on a high contrast feature of the object you are imaging. Focusing on small features such as a moon’s shadow on Jupiter or Cassini’s division in the ring of Saturn will guarantee best focus across the entire image.Once the frames of your video are stacked, the overall brightness of the composite (stacked) image is usually brighter than its individual component frames. For this reason it is best to keep the brightness of the video image seen on the screen dimmer than you would normally desire. It is important thatno part of the image is over-exposed to assure the maximum amount of detail in the final composite image.CollimationNo matter what type of telescope you image with, poor collimation (alignment of the optics) willruin your chances for a good image. Before you begin imaging, always check the collimation of your instrument and make adjustments if necessary. Refer to your telescopes owner’s manual for instruction in collimating the optics.Finding ObjectsAt first it can be difficult to locate individual planets due to their relative brightness. To make it easier to initially find your object in the imaging window, increase the brightness and gain controls onthe Exposure Tool Bar. This will allow you to better see the object as it passes through the imaging window. Once the object is located and centered, you can adjust the setting until the object is at the desired brightness and contrast.how long to take streaming VideoAt first you may think that the more frames you record the better. However there are some limitation to the duration of video and the amount of frames you can acquire. Resolution and file size can both limit the length of time of your video.Since the NexImage will combine as many sharp frames as possible to achieve one high quality image, you don’t want to take so many images that you start to detect the rotation of the planet, especially on Jupiter which makes one complete rotation in under 10 hours!Also each frame of a high resolution video can equal a large file size. Hundreds of frames cantake up much of your hard drive space. Since file sizes of each video taken can be quite large, it is recommended that you save your video data onto a CD-ROM. This way you can have a library of files stored for processing without filling up your hard drive in the process.9AutoguidingYour NexImage 5 can be used as an autoguiding camera when used in junction with an off-axis guider or piggyback guidescope. NexImage is compatible with many autoguiding software such as MetaGuide (/Bliss/MetaGuide), and a guider port interface such as GPUSB from Shoestring Astronomy ().Specifications:Sensor 1 /2.5” format, color CMOS Sensor (Micron MT9P031)Camera Resolution5MP (2592 x 1944)Sensor Size 5.7mm x 4.28mm (7 mm diagonal)Pixel Size 2.2 micron squareSensitivity 1.4 V/ lux-sec (@550nm)USB Cable High-Speed USB 2.0 cable102835 Columbia Street • Torrance, CA 90503 U.S.A.Telephone: 310.328.9560 • Fax: 310.212.5835©2012 CelestronAll rights reserved. • Printed in China • 03-12FCC Note: This equipment has been tested and found to comply with the limits for a Class B digital device, pursuant to part 15 of the FCC Rules. These limits are designed to provide reasonable protection against harmful interference in a residential installation. This equipment generates, uses, and can radiate radio frequency energy and, if not installed and used in accordance with the instructions, may cause harmful interference to radio communications. However, there is no guarantee that interference will not occur in a particular installation. If this equipment does cause harmful interference to radio or television reception, which can be determined by turning the equipment off and on, the user is encouraged to try to correct the interference by one or more of the following measures:• Reorient or relocate the receiving antenna.• Increase the separation between the equipment and receiver.• Connect the equipment into an outlet on a circuit different from that to which the receiver is connected.• Consult the dealer or an experienced radio/TV technician for help.Designed and intended for those 13 years of age and older.。
Method of Face Recognition Based on Red-BlackWavelet Transform and PCAYuqing He, Huan He, and Hongying YangDepartment of Opto-Electronic Engineering,Beijing Institute of Technology, Beijing, P.R. China, 10008120701170@。
cnAbstract。
With the development of the man—machine interface and the recogni—tion technology, face recognition has became one of the most important research aspects in the biological features recognition domain. Nowadays, PCA(Principal Components Analysis) has applied in recognition based on many face database and achieved good results. However, PCA has its limitations: the large volume of computing and the low distinction ability。
In view of these limitations, this paper puts forward a face recognition method based on red—black wavelet transform and PCA. The improved histogram equalization is used to realize image pre-processing in order to compensate the illumination. Then, appling the red—black wavelet sub—band which contains the information of the original image to extract the feature and do matching。
904N.E.Huang and others10.Discussion98711.Conclusions991References993 A new method for analysing has been devel-oped.The key part of the methodany complicated data set can be decomposed intoof‘intrinsic mode functions’Hilbert trans-This decomposition method is adaptive,and,highly efficient.Sinceapplicable to nonlinear and non-stationary processes.With the Hilbert transform,Examplesthe classical nonlinear equation systems and dataare given to demonstrate the power new method.data are especially interesting,for serve to illustrate the roles thenonlinear and non-stationary effects in the energy–frequency–time distribution.Keywords:non-stationary time series;nonlinear differential equations;frequency–time spectrum;Hilbert spectral analysis;intrinsic time scale;empirical mode decomposition1.Introductionsensed by us;data analysis serves two purposes:determine the parameters needed to construct the necessary model,and to confirm the model we constructed to represent the phe-nomenon.Unfortunately,the data,whether from physical measurements or numerical modelling,most likely will have one or more of the following problems:(a)the total data span is too short;(b)the data are non-stationary;and(c)the data represent nonlinear processes.Although each of the above problems can be real by itself,the first two are related,for a data section shorter than the longest time scale of a sta-tionary process can appear to be non-stationary.Facing such data,we have limited options to use in the analysis.Historically,Fourier spectral analysis has provided a general method for examin-the data analysis has been applied to all kinds of data.Although the Fourier transform is valid under extremely general conditions(see,for example,Titchmarsh1948),there are some crucial restrictions of Proc.R.Soc.Lond.A(1998)Nonlinear and non-stationary time series analysis905the Fourier spectral analysis:the system must be linear;and the data must be strict-ly periodic or stationary;otherwise,the resulting spectrum will make little physicalsense.to the Fourier spectral analysis methods.Therefore,behoves us review the definitions of stationarity here.According to the traditional definition,a time series,X (t ),is stationary in the wide sense,if,for all t ,E (|X (t )2|)<∞,E (X (t))=m,C (X (t 1),X (t 2))=C (X (t 1+τ),X (t 2+τ))=C (t 1−t 2),(1.1)in whichE (·)is the expected value defined as the ensemble average of the quantity,and C (·)is the covariance function.Stationarity in the wide sense is also known as weak stationarity,covariance stationarity or second-order stationarity (see,forexample,Brockwell &Davis 1991).A time series,X (t ),is strictly stationary,if the joint distribution of [X (t 1),X (t 2),...,X (t n )]and [X (t 1+τ),X (t 2+τ),...,X (t n +τ)](1.2)are the same for all t i and τ.Thus,a strictly stationaryprocess with finite second moments is alsoweakly stationary,but the inverse is not true.Both definitions arerigorous but idealized.Other less rigorous definitions have also beenused;for example,that is stationary within a limited timespan,asymptotically stationary is for any random variableis stationary when τin equations (1.1)or (1.2)approaches infinity.In practice,we can only have data for finite time spans;these defini-tions,we haveto makeapproximations.Few of the data sets,from either natural phenomena or artificial sources,can satisfy these definitions.It may be argued thatthe difficulty of invoking stationarity as well as ergodicity is not on principlebut on practicality:we just cannot have enough data to cover all possible points in thephase plane;therefore,most of the cases facing us are transient in nature.This is the reality;we are forced to face it.Fourier spectral analysis also requires linearity.can be approximated by linear systems,the tendency tobe nonlinear whenever their variations become finite Compounding these complications is the imperfection of or numerical schemes;theinteractionsof the imperfect probes even with a perfect linear systemcan make the final data nonlinear.For the above the available data are ally of finite duration,non-stationary and from systems that are frequently nonlinear,either intrinsicallyor through interactions with the imperfect probes or numerical schemes.Under these conditions,Fourier spectral analysis is of limited use.For lack of alternatives,however,Fourier spectral analysis is still used to process such data.The uncritical use of Fourier spectral analysis the insouciant adoption of the stationary and linear assumptions may give cy range.a delta function will giveProc.R.Soc.Lond.A (1998)906N.E.Huang and othersa phase-locked wide white Fourier spectrum.Here,added to the data in the time domain,Constrained bythese spurious harmonics the wide frequency spectrum cannot faithfully represent the true energy density in the frequency space.More seri-ously,the Fourier representation also requires the existence of negative light intensity so that the components can cancel out one another to give thefinal delta function. Thus,the Fourier components might make mathematical sense,but do not really make physical sense at all.Although no physical process can be represented exactly by a delta function,some data such as the near-field strong earthquake records areFourier spectra.Second,tions;wave-profiles.Such deformations,later,are the direct consequence of nonlinear effects.Whenever the form of the data deviates from a pure sine or cosine function,the Fourier spectrum will contain harmonics.As explained above, both non-stationarity and nonlinearity can induce spurious harmonic components that cause energy spreading.The consequence is the misleading energy–frequency distribution forIn this paper,modemode functions The decomposition is based on the direct extraction of theevent on the time the frequency The decomposition be viewed as an expansion of the data in terms of the IMFs.Then,based on and derived from the data,can serve as the basis of that expansion linear or nonlinear as dictated by the data,Most important of all,it is adaptive.As will locality and adaptivity are the necessary conditions for the basis for expanding nonlinear and non-stationary time orthogonality is not a necessary criterionselection for a nonlinearon the physical time scaleslocal energy and the instantaneous frequencyHilbert transform can give us a full energy–frequency–time distribution of the data. Such a representation is designated as the Hilbert spectrum;it would be ideal for nonlinear and non-stationary data analysis.We have obtained good results and new insights by applying the combination of the EMD and Hilbert spectral analysis methods to various data:from the numerical results of the classical nonlinear equation systems to data representing natural phe-nomena.The classical nonlinear systems serve to illustrate the roles played by the nonlinear effects in the energy–frequency–time distribution.With the low degrees of freedom,they can train our eyes for more complicated cases.Some limitations of this method will also be discussed and the conclusions presented.Before introducing the new method,we willfirst review the present available data analysis methods for non-stationary processes.Proc.R.Soc.Lond.A(1998)Nonlinear and non-stationary time series analysis9072.Review of non-stationary data processing methodsWe willfirstgivea brief survey of themethodsstationary data.are limited to linear systems any method is almost strictly determined according to the special field in which the application is made.The available methods are reviewed as follows.(a )The spectrogramnothing but a limited time window-width Fourier spectral analysis.the a distribution.Since it relies on the tradition-al Fourier spectral analysis,one has to assume the data to be piecewise stationary.This assumption is not always justified in non-stationary data.Even if the data are piecewise stationary how can we guarantee that the window size adopted always coincides with the stationary time scales?What can we learn about the variations longer than the local stationary time scale?Will the collection of the locally station-ary pieces constitute some longer period phenomena?Furthermore,there are also practical difficulties in applying the method:in order to localize an event in time,the window width must be narrow,but,on the other hand,the frequency resolu-tion requires longer time series.These conflicting requirements render this method of limited usage.It is,however,extremely easy to implement with the fast Fourier transform;thus,ithas attracted a wide following.Most applications of this methodare for qualitative display of speech pattern analysis (see,for example,Oppenheim &Schafer 1989).(b )The wavelet analysisThe wavelet approach is essentially an adjustable window Fourier spectral analysiswith the following general definition:W (a,b ;X,ψ)=|a |−1/2∞−∞X (t )ψ∗ t −b ad t,(2.1)in whichψ∗(·)is the basic wavelet function that satisfies certain very general condi-tions,a is the dilation factor and b is the translationof theorigin.Although time andfrequency do not appear explicitly in the transformed result,the variable 1/a givesthe frequency scale and b ,the temporal location of an event.An intuitive physical explanation of equation (2.1)is very simple:W (a,b ;X,ψ)is the ‘energy’of X ofscale a at t =b .Because of this basic form of at +b involvedin thetransformation,it is also knownas affinewavelet analysis.For specific applications,the basic wavelet function,ψ∗(·),can be modified according to special needs,but the form has to be given before the analysis.In most common applications,however,the Morlet wavelet is defined as Gaussian enveloped sine and cosine wave groups with 5.5waves (see,for example,Chan 1995).Generally,ψ∗(·)is not orthogonalfordifferent a for continuous wavelets.Although one can make the wavelet orthogonal by selecting a discrete set of a ,thisdiscrete wavelet analysis will miss physical signals having scale different from theselected discrete set of a .Continuous or discrete,the wavelet analysis is basically a linear analysis.A very appealing feature of the wavelet analysis is that it provides aProc.R.Soc.Lond.A (1998)908N.E.Huang and othersuniform resolution for all the scales.Limited by the size of thebasic wavelet function,the downside of the uniform resolution is uniformly poor resolution.Although wavelet analysis has been available only in the last ten years or so,it hasbecome extremelypopular.Indeed,it is very useful in analysing data with gradualfrequency changes.Since it has an analytic form for the result,it has attracted extensive attention of the applied mathematicians.Most of its applications have been in edge detection and image compression.Limited applications have also been made to the time–frequency distribution in time series (see,for example,Farge 1992;Long et al .1993)andtwo-dimensionalimages (Spedding et al .1993).Versatile as the wavelet analysis is,the problem with the most commonly usedMorlet wavelet is its leakage generated by the limited length of the basic wavelet function,whichmakesthe quantitativedefinitionof the energy–frequency–time dis-tribution difficult.Sometimes,the interpretation of the wavelet can also be counter-intuitive.For example,to define a change occurring locally,one must look for theresult in the high-frequencyrange,for the higher the frequency the more localized thebasic wavelet will be.If a local event occurs only in the low-frequency range,one willstill be forced to look for its effects inthe high-frequencyrange.Such interpretationwill be difficultif it is possible at all (see,for example,Huang et al .1996).Another difficulty of the wavelet analysis is its non-adaptive nature.Once the basic waveletis selected,one will have to use it to analyse all the data.Since the most commonlyused Morlet wavelet is Fourier based,it also suffers the many shortcomings of Fouri-er spectral analysis:it can only give a physically meaningful interpretation to linear phenomena;it can resolve the interwave frequency modulation provided the frequen-cy variationis gradual,but it cannot resolve the intrawave frequency modulation because the basic wavelet has a length of 5.5waves.Inspite of all these problems,wavelet analysisisstillthe bestavailable non-stationary data analysis method so far;therefore,we will use it in this paper as a reference to establish the validity and thecalibration of the Hilbert spectrum.(c )The Wigner–Ville distributionThe Wigner–Ville distribution is sometimes alsoreferred toas the Heisenberg wavelet.By definition,it is the Fourier transform of the central covariance function.For any time series,X (t ),we can define the central variance as C c (τ,t )=X (t −12τ)X ∗(t +12τ).(2.2)Then the Wigner–Ville distribution is V (ω,t )=∞−∞C c (τ,t )e −i ωτd τ.(2.3)This transform has been treated extensively by Claasen &Mecklenbr¨a uker (1980a ,b,c )and by Cohen (1995).It has been extremely popular with the electrical engi-neering community.The difficulty with this method is the severe cross terms as indicated by the exis-tence of negativepowerfor some frequency ranges.Although this shortcoming canbe eliminated by using the Kernel method (see,for example,Cohen 1995),the resultis,then,basically that of a windowed Fourier analysis;therefore,itsuffers all thelim-itations of the Fourier analysis.An extension of this method has been made by Yen(1994),who used the Wigner–Ville distribution to define wave packets that reduce Proc.R.Soc.Lond.A (1998)Nonlinear and non-stationary time series analysis909 a complicated data set to afinite number of simple components.This extension is very powerful and can be applied to a variety of problems.The applications to complicated data,however,require a great amount of judgement.(d)Evolutionary spectrumThe evolutionary spectrum wasfirst proposed by Priestley(1965).The basic idea is to extend the classic Fourier spectral analysis to a more generalized basis:from sine or cosine to a family of orthogonal functions{φ(ω,t)}indexed by time,t,and defined for all realω,the frequency.Then,any real random variable,X(t),can beexpressed asX(t)= ∞−∞φ(ω,t)d A(ω,t),(2.4)in which d A(ω,t),the Stieltjes function for the amplitude,is related to the spectrum asE(|d A(ω,t)|2)=dµ(ω,t)=S(ω,t)dω,(2.5) whereµ(ω,t)is the spectrum,and S(ω,t)is the spectral density at a specific time t,also designated as the evolutionary spectrum.If for eachfixedω,φ(ω,t)has a Fourier transformφ(ω,t)=a(ω,t)e iΩ(ω)t,(2.6) then the function of a(ω,t)is the envelope ofφ(ω,t),andΩ(ω)is the frequency.If, further,we can treatΩ(ω)as a single valued function ofω,thenφ(ω,t)=α(ω,t)e iωt.(2.7) Thus,the original data can be expanded in a family of amplitude modulated trigono-metric functions.The evolutionary spectral analysis is very popular in the earthquake communi-ty(see,for example,Liu1970,1971,1973;Lin&Cai1995).The difficulty of its application is tofind a method to define the basis,{φ(ω,t)}.In principle,for this method to work,the basis has to be defined a posteriori.So far,no systematic way has been offered;therefore,constructing an evolutionary spectrum from the given data is impossible.As a result,in the earthquake community,the applications of this method have changed the problem from data analysis to data simulation:an evo-lutionary spectrum will be assumed,then the signal will be reconstituted based on the assumed spectrum.Although there is some general resemblance to the simulated earthquake signal with the real data,it is not the data that generated the spectrum. Consequently,evolutionary spectrum analysis has never been very useful.As will be shown,the EMD can replace the evolutionary spectrum with a truly adaptive representation for the non-stationary processes.(e)The empirical orthogonal function expansion(EOF)The empirical orthogonal function expansion(EOF)is also known as the principal component analysis,or singular value decomposition method.The essence of EOF is briefly summarized as follows:for any real z(x,t),the EOF will reduce it toz(x,t)=n1a k(t)f k(x),(2.8)Proc.R.Soc.Lond.A(1998)910N.E.Huang and othersin whichf j·f k=δjk.(2.9)The orthonormal basis,{f k},is the collection of the empirical eigenfunctions defined byC·f k=λk f k,(2.10)where C is the sum of the inner products of the variable.EOF represents a radical departure from all the above methods,for the expansion basis is derived from the data;therefore,it is a posteriori,and highly efficient.The criticalflaw of EOF is that it only gives a distribution of the variance in the modes defined by{f k},but this distribution by itself does not suggest scales or frequency content of the signal.Although it is tempting to interpret each mode as indepen-dent variations,this interpretation should be viewed with great care,for the EOF decomposition is not unique.A single component out of a non-unique decomposition, even if the basis is orthogonal,does not usually contain physical meaning.Recently, Vautard&Ghil(1989)proposed the singular spectral analysis method,which is the Fourier transform of the EOF.Here again,we have to be sure that each EOF com-ponent is stationary,otherwise the Fourier spectral analysis will make little sense on the EOF components.Unfortunately,there is no guarantee that EOF compo-nents from a nonlinear and non-stationary data set will all be linear and stationary. Consequently,singular spectral analysis is not a real improvement.Because of its adaptive nature,however,the EOF method has been very popular,especially in the oceanography and meteorology communities(see,for example,Simpson1991).(f)Other miscellaneous methodsOther than the above methods,there are also some miscellaneous methods such as least square estimation of the trend,smoothing by moving averaging,and differencing to generate stationary data.Methods like these,though useful,are too specialized to be of general use.They will not be discussed any further here.Additional details can be found in many standard data processing books(see,for example,Brockwell &Davis1991).All the above methods are designed to modify the global representation of the Fourier analysis,but they all failed in one way or the other.Having reviewed the methods,we can summarize the necessary conditions for the basis to represent a nonlinear and non-stationary time series:(a)complete;(b)orthogonal;(c)local;and (d)adaptive.Thefirst condition guarantees the degree of precision of the expansion;the second condition guarantees positivity of energy and avoids leakage.They are the standard requirements for all the linear expansion methods.For nonlinear expansions,the orthogonality condition needs to be modified.The details will be discussed later.But even these basic conditions are not satisfied by some of the above mentioned meth-ods.The additional conditions are particular to the nonlinear and non-stationary data.The requirement for locality is the most crucial for non-stationarity,for in such data there is no time scale;therefore,all events have to be identified by the time of their occurences.Consequently,we require both the amplitude(or energy) and the frequency to be functions of time.The requirement for adaptivity is also crucial for both nonlinear and non-stationary data,for only by adapting to the local variations of the data can the decomposition fully account for the underlying physics Proc.R.Soc.Lond.A(1998)Nonlinear and non-stationary time series analysis911of the processes and not just to fulfil the mathematical requirements for fitting the data.This is especially important for the nonlinear phenomena,for a manifestation of nonlinearity is the ‘harmonic distortion’in the Fourier analysis.The degree of distortion depends on the severity of nonlinearity;therefore,one cannot expect a predetermined basis to fit all the phenomena.An easy way to generate the necessary adaptive basis is to derive the basis from the data.In this paper,we will introduce a general method which requires two steps in analysing the data as follows.The first step is to preprocess the data by the empirical mode decomposition method,with which the data are decomposed into a number of intrinsic mode function components.Thus,we will expand the data in a basis derived from the data.The second step is to apply the Hilbert transform to the decomposed IMFs and construct the energy–frequency–time distribution,designated as the Hilbert spectrum,from which the time localities of events will be preserved.In other words,weneed the instantaneous frequency and energy rather than the global frequency and energy defined by the Fourier spectral analysis.Therefore,before goingany further,we have to clarify the definition of the instantaneous frequency.3.Instantaneous frequencyis to accepting it only for special ‘monocomponent’signals 1992;Cohen 1995).Thereare two basicdifficulties with accepting the idea of an instantaneous fre-quency as follows.The first one arises from the influence of theFourier spectral analysis.In the traditional Fourier analysis,the frequency is defined for thesineor cosine function spanning the whole data length with constant ampli-tude.As an extension of this definition,the instantaneous frequencies also have torelate to either a sine or a cosine function.Thus,we need at least one full oscillationof a sineor a cosine wave to define the local frequency value.According to this logic,nothing full wave will do.Such a definition would not make sense forThe secondarises from the non-unique way in defining the instantaneousfrequency.Nevertheless,this difficulty is no longer serious since the introduction ofthe meanstomakethedata analyticalthrough the Hilbert transform.Difficulties,however,still exist as ‘paradoxes’discussed by Cohen (1995).For an arbitrary timeseries,X (t ),we can always have its Hilbert Transform,Y (t ),as Y (t )=1πP∞−∞X (t )t −t d t,(3.1)where P indicates the Cauchy principal value.This transformexists forallfunctionsof class L p(see,for example,Titchmarsh 1948).With this definition,X (t )and Y (t )form the complex conjugate pair,so we can have an analytic signal,Z (t ),as Z (t )=X (t )+i Y (t )=a (t )e i θ(t ),(3.2)in which a (t )=[X 2(t )+Y 2(t )]1/2,θ(t )=arctanY (t )X (t ).(3.3)Proc.R.Soc.Lond.A (1998)912N.E.Huang andothers Theoretically,there are infinitely many ways of defining the imaginary part,but the Hilbert transform provides a unique way of defining the imaginary part so that the result is ananalyticfunction.A brief tutorial on the Hilbert transform with theemphasis on its physical interpretation can be found in Bendat &Piersol is the bestlocal fitan amplitude and phase varying trigonometric function to X (t ).Even with the Hilbert transform,there is still controversy in defining the instantaneous frequency as ω=d θ(t )d t .(3.4)This leads Cohen (1995)to introduce the term,‘monocomponent function’.In prin-ciple,some limitations on the data are necessary,forthe instantaneous frequencygiven in equation (3.4)is a single value function of time.At any given time,thereis only one frequency value;therefore,it can only represent one component,hence ‘monocomponent’.Unfortunately,no cleardefinition of the ‘monocomponent’signalwas given to judge whether a function is or is not ‘monocomponent’.For lack ofa precise definition,‘narrow band’was adopted a on the data for the instantaneous frequency to make sense (Schwartz et al .1966).There are two definitions for bandwidth.The first one is used in the study of the probability properties of the signalsand waves,wherethe processes are assumed tobe stationary and Gaussian.Then,the bandwidth can be defined in spectral moments The expected number of zero crossings per unit time is given byN 0=1π m 2m 0 1/2,(3.5)while the expected number of extrema per unit time is given byN 1=1π m 4m 2 1/2,(3.6)in which m i is the i th moment of the spectrum.Therefore,the parameter,ν,definedas N 21−N 20=1π2m 4m 0−m 22m 2m 0=1π2ν2,(3.7)offers a standard bandwidth measure (see,for example,Rice 1944a,b ,1945a,b ;Longuet-Higgins 1957).For a narrow band signal ν=0,the expected numbers extrema and zero crossings have to equal.the spectrum,but in a different way.coordinates as z (t )=a (t )e i θ(t ),(3.8)with both a (t )and θ(t )being functions of time.If this function has a spectrum,S (ω),then the mean frequency is given byω = ω|S (ω)|2d ω,(3.9)Proc.R.Soc.Lond.A (1998)Nonlinear and non-stationary time series analysis913which can be expressed in another way asω =z ∗(t )1i dd tz (t )d t=˙θ(t )−i ˙a (t )a (t )a 2(t )d t =˙θ(t )a 2(t )d t.(3.10)Based on this expression,Cohen (1995)suggested that ˙θbe treated as the instanta-neous frequency.With these notations,the bandwidth can be defined asν2=(ω− ω )2 ω 2=1 ω 2(ω− ω )2|S (ω)|2d ω=1 ω 2z ∗(t ) 1i d d t− ω 2z (t )d t =1 ω 2 ˙a 2(t )d t +(˙θ(t )− ω )2a 2(t )d t .(3.11)For a narrow band signal,this value has to be small,then both a and θhave to begradually varying functions.Unfortunately,both equations (3.7)and (3.11)defined the bandwidth in the global sense;they are both overly restrictive and lack preci-sion at the same time.Consequently,the bandwidth limitation on the Hilbert trans-form to give a meaningful instantaneous frequency has never been firmly established.For example,Melville (1983)had faithfully filtered the data within the bandwidth requirement,but he still obtained many non-physical negative frequency values.It should be mentioned here that using filtering to obtain a narrow band signal is unsat-isfactory for another reason:the filtered data have already been contaminated by the spurious harmonics caused by the nonlinearity and non-stationarity as discussed in the introduction.In order to obtain meaningful instantaneous frequency,restrictive conditions have to be imposed on the data as discussed by Gabor (1946),Bedrosian (1963)and,more recently,Boashash (1992):for any function to have a meaningful instantaneous frequency,the real part of its Fourier transform has to have only positive frequency.This restriction can be proven mathematically as shown in Titchmarsh (1948)but it is still global.For data analysis,we have to translate this requirement into physically implementable steps to develop a simple method for applications.For this purpose,we have to modify the restriction condition from a global one to a local one,and the basis has to satisfy the necessary conditions listed in the last section.Let us consider some simple examples to illustrate these restrictions physically,by examining the function,x (t )=sin t.(3.12)Its Hilbert transform is simply cos t .The phase plot of x –y is a simple circle of unit radius as in figure 1a .The phase function is a straight line as shown in figure 1b and the instantaneous frequency,shown in figure 1c ,is a constant as expected.If we move the mean offby an amount α,say,then,x (t )=α+sin t.(3.13)Proc.R.Soc.Lond.A (1998)。
计算技术与自动化Computing Technology and Automation第40卷第1期2 0 2 1年3月Vol. 40,No. 1Mar. 2 02 1文章编号:1003-6199( 2021 )01-0101 — 03DOI : 10. 16339/j. cnki. jsjsyzdh. 202101019基于3次B 样条小波变换的改进自适应阈值边缘检测算法王 煜J 谢 政,朱淳钊,夏建高(湖北工程职业学院建筑与环境艺术学院,湖北黄石435005)摘要:针对含噪声图像边缘提取问题,提出了一种改进NormalShrink 自适应阈值去噪算法。
该算法首先通过小波变换和局部模极大值法提取出可能包含图像边缘特征的小波系数,利用边缘像素之间特殊的空间关系以及噪声在各级小波分解尺度下的不同效应,构建适合各个尺度级的改进NormalShrink 自适应阈值,并依此对提取出的小波系数进行筛选。
实验结果表明,与改进的Candy 算子和传统的NormalShrink 自 适应阈值相比,本方法提取出的图像边缘较为完整清晰,峰值信噪比提升约6 db o关键词:边缘提取;小波变换;自适应阈值;峰值信噪比中图分类号:TP312文献标识码:AAn Improved Adaptive Threshold Edge Detection AlgorithmBased on Cubic B-spline Wavelet TransformWANG Yu f , XIE Zheng,ZHU Chun-zhao ,XIA Jian-gao(School of Architecture and Environmental Art, Hubei Engineering Institute, Huangshi, Hubei 435005, China)Abstract : In order to solve the problem of noisy image edge detection, an improved NormalShrink adaptive waveletthreshold is put forward on the foundation of combining edge detection and denoising . According to the different characteris tics of noise at different wavelet scales and the special spatial relationship between the edge pixels , the algorithm first extract wavelet coefficients which may contain image edge feature by using wavelet transform and local maximum mode, and thenconstruct an improved NormalShrink adaptive threshold of each scale level which is used to select the extracted wavelet coef ficients. Experimental results show that this method can keep imagers edges clear and increase PSNR about 6 db.Key words :edge detection ; wavelet transform ; adaptive threshold ; PSNR图像边缘信息的识别和提取在图像分割、图像 识别等领域有着重要的应用,提取出清晰有效的边缘是一个热点研究方向。
引用格式:苏莹莹, 毛海旭. 小波变换和CNN 涡旋压缩机故障诊断[J]. 中国测试,2023, 49(4): 92-97. SU Yingying, MAO Haixu.Fault diagnosis of scroll compressor based on wavelet transform and CNN[J]. China Measurement & Test, 2023, 49(4): 92-97. DOI :10.11857/j.issn.1674-5124.2021070062小波变换和CNN 涡旋压缩机故障诊断苏莹莹, 毛海旭(沈阳大学机械工程学院,辽宁 沈阳 110000)摘 要: 针对传统单尺度信号分析难以有效解决涡旋压缩机故障诊断中的故障特征信息多尺度耦合问题,提出一种基于小波变换和卷积神经网络的涡旋压缩机故障诊断方法。
首先将采集到的振动信号进行连续小波变换生成时频图,并对时频图进行网格化规范处理,将预处理后的时频图作为特征图输入Alexnet 卷积神经网络,通过不断调节网络参数,得出最为理想的神经网络模型,以此实现对涡旋压缩机故障类型的辨识诊断。
结果表明,该方法针对涡旋压缩机故障类型的识别准确率达到94.6%,与传统多尺度排列熵、信息熵熵距的故障诊断方法相比,该故障识别方法具有更高的准确率。
关键词: 故障诊断; 振动信号; 小波变换; 卷积神经网络中图分类号: U226.8+1;TB9文献标志码: A文章编号: 1674–5124(2023)04–0092–06Fault diagnosis of scroll compressor based on wavelet transform and CNNSU Yingying, MAO Haixu(School of Mechanical Engineering, Shenyang University, Shenyang 110000, China)Abstract : In order to solve the problem that traditional single-scale signal analysis is difficult to effectively solve the problem of multi-scale coupling of fault feature information in the fault diagnosis of scroll compressors, a fault diagnosis method based on wavelet transform and convolutional neural network(CNN) is proposed. Firstly, the collected vibration signal is analyzed by continuous wavelet transform to generate time-frequency diagram. And the generated time-frequency diagram is gridded and normalized. Then, it has to be inputtd to Alexnet convolutional neural network, and the network parameters are adjusted to obtain the most ideal network model, so as to realize the identification and diagnosis fault types of scroll compressors. The results show that the recognition accuracy reaches 94.6%, with higher accuracy than the traditional methods of multi-scale permutation entropy and distance of the information entropy.Keywords : fault diagnosis; vibration signal; wavelet transform; convolutional neural networks0 引 言涡旋压缩机因为其体积小、集成度高、运行稳定、噪声低等优点,在工业、农业、交通运输等众多领域广泛应用。
武汉大学硕士学位论文基于小波变换的周期信号的特征提取姓名:***申请学位级别:硕士专业:应用数学指导教师:***20050401摘要时颓分析旨在构造一种时闯稆焱率翡密度函数,激搦示信号中所包含静频率分楚及演化特性。
在时频分析方法中,小波交换是近年来迅速发展起来的一种新的正具,它在时域和频域上同时县有良好的局部化特性,并可产生自适应的时频窗,无论分析低频或高频局部信号,它都能自动调节时频寓,以适应实际分析的需鼹。
小波变换程时频分横中具有很强的烫活性,从而可以聚焦到信号对段驷频段的任懑细节,被人们蔻炎时频分糖浆显微镶。
程瓣频分撰中一般要获取信号豹特薤倍感,逶j霪对绩号俸交换后,信号鼢将征信意搿瑷懋可煞的突显出来敬稻于特征提取。
由于信弩的奇异点能够翔画信号的细节,通过Lipschitz指数来袭征信号的奇异性,然后根据模极大值与Lipschitz指数的关系,就可以通过模极大值来提取出这毖奇异点,这就是基于小波变换的模极大饿的特征提取。
为此,本文首先阐述了时频分柝技术及研究现状,并燎现在鬻用魄时频分板方法进行了分类,对不同弱对频分橱技术撵了麓要奔缨,著对其键缺点、彼我阂静辖互关系滋孬了较为详赣豹论述。
然后溺述了,j、波交换靛基本鹫论和模校大值原疆。
针对当前普遍采餍韵基于二进小波变换瀚模极大值特征提敬方法的不足,建立了蒸于连续小波变换的模极大值特征提取方法,具有更细致的局部模极大值刻画能力。
最后本文通过确定周期信号连续小波变换的模极大值和极点的位鼹来定魑的描述信号的基本信息,如振幅、频率、棚位、瞬时频率和瞬时振螭等。
关键词:特征提取Lipschitz指数模极大值连续小波变换ABSTRACTTheaimoftime-丘equenayanalysisistofindadensityintimeandfrequencythatindicateswhichfrequenciesarepresentinthesignalandhowtheyevolvewithtime.Inthegeneralmethodoftime一丘equencyanalysis,waveletanalysisisintemationalyrccongnizeduptotheminutetoolsforanalysingtime-frequenc孓Itischieflyduetothe“adaptivefeature’’and“mathematicalmierotelescopefeature”.waveletanalysisisbecomingafocuspointofmanysciences,andisfondlydelightedastoolsforSOman),scientificworkers.Itplaysanimportantroleinthesignal&informationprocessing.Thepointofsharpvariationorsingularitiesusuallycarrytheimportantinformationaboutthesignals.Theadventofwavelettransformprovidesasuitableframeworkforstudyingthemultiscaletransientrepresentationofsignals.Inmathematics,singularityCallbedescribedbylipschitzexponents.MallathasprovedthatthelipschitzexponentsCanbeobtainedbycalculatingthedecayofthelocalmodulemaximumofthewavelettransform,whichiSthemodulemaximumfeatureextractionmethodsbasedonwavelettransform.Thispapersummarizesbrieflydifferenttime-frquencyanalysistechniquesandresearchoftime·frequencyanalysis.Itdiscussesindetailtheadvantagesanddisadvantagesofthesetechniquesandtheirrelations。
6.小波分析 wavelet analysis自从学习过佛瑞艾尔变形和频率估计,作者对小波分析产生的兴趣,开始阅读一些相关的资料。
同时发现正在给自己上数学课的David Walnut 教授是小波分析的泰山北斗。
他写的这方面的教课书遍布全球。
作者的心中又产生的无比的崇敬与羡慕,所以将当前世界最先进的小波分析技巧写出来与大家分享。
读过《R语言时间系列中文教程》都应该知道如何使用弗瑞艾尔变形估计频率。
但必须假设,被估计的频率是始终存在于波动之中的。
更经常的状况是在一整个波中某一频率只在这个波中的一小部分出现。
使用弗瑞艾尔变形不可能监测到在某一时间点上的频率变化,因为它假设所估计的频率都是自始至终存在的。
例如,下面的波中有一个很慢的频率是始终存在的,在中间部分突然出现了频率非常高的新波动,而且很快就消失了。
这样的波动需要使用小波分析。
t=1:500c1 = 2*cos(2*pi*t/150 + .6*pi)plot.ts(c1)t2= ifelse(t>200 & t<300,t,0)c2=1*cos(2*pi*t2/10 + .6*pi)par(new=T)plot.ts(c2)par(new=F)例如心电图是有来测量人心脏跳动的手段,心脏的各个组成部分不是同时不停的工作的,在一个心房收缩的过程中某些的心房是休息的。
所以在心电图上来看,微小的波动是突然出现仅仅持续很短的时间就消失了,像这样的微小波动就要使用小波分析来捕捉。
心电图是医生进行临床诊断的方法之一,通过阅读心电图医生可以推测病人的心脏是正常的还是哪里出现的毛病。
接下来我们要介绍使用小波分析计算机智能进行自动诊断,也就是说我们积累了很多人的心电图数据通过小波分析将心电图的特色提取出来。
这些特色将告诉我们病人的心脏是正常的还是有状况。
通过这些心电图的数据我们可以建立一个数学模型,当有一个新的病人进来我们就可以对他的心电图进行诊断。
Signal filtering, Wavelet transform, Hilbert-Huang transformEn jämförelseNazar Dino mohammedSammanfattningDetta examensarbete behandlar och jämför DWT (Discrete Wavelet Transform) och HHT(Hilbert-Huang Transform). I många sammanhangarer och tillämpningar, tid- frekvens representation av en signal är viktig och betydelsefull. Tid- frekvens analyser ger de mest betydelsefulla informationer om signalen frekvens, deras styrka och hur de ändras sig med tiden.Traditionella metoder som Fouier transform, kan ge bara en frekvens baserade representation av signalen och kan tillämpas bara på linjära processer. De flesta processer i verkligheten är icke linjär och icke stationära processer och därför andra metoder som kan tillämpas på alla signaler och ger en tid- frekvens approximation av signaler, är nödvändiga. DWT och HHT är några av de metoder som kan användas för detta syfte.Ett annat viktig del av signal processning handlar om klassificering av signaler ochkänna igen signaler och detta kan göras med PCA (Principal component analysis) och NN (Neural network). Men innan dess måste viktigaste information om tid och frekvens av signaler bestämmas och vi kan göra detta med DWT och HHT. Vi tarnågra riktiga signaler (infrasound signaler)och tillämpar de olika metoder och jämför dem.AbstractIn this work, we will compare the discrete wavelet transform (DWT) and Hilbert Huang transform (HHT) as pre-processing techniques for signal classification using Neural Network. Time-frequency analysis is important part of signal processing. Time-frequency analysis is the process of determining what frequencies are present in a signal, how strong they are and how they change over time. DWT and HHT techniques, determine the different frequencies in the signal.Wavelet transform is common approach to study time-frequency resolution in signal processing. Wavelet transform can be applied to linear and non stationary processes while HHT can be applied to non linear processes and non stationary process.The physical processes in the nature are usually non-linear and non-stationary.For example infrasound signals can be considered non-linear and non-stationary processes. In our experiments, we use infrasound signals from different events. Classification and recognition of signals are important parts of signal processing and they have many application such face recognition, fingerprint matching, voice recognition, to name a few. This can be done for example by methods like Neural networks (NN), however before that, we need to take the most significant information from signals like time and frequency. Both DWT and HHT are optimal methods for this purpose .In this work we use some real signals (infrasound signals) and processing them with DWT, HHT and NN.1 IntroductionThe goal of our project is to compare the discrete wavelet transform (DWT) and Hilbert Huang transform (HHT) as pre-processing techniques for signal classification using Neural Network. We use infrasound signals from two events: sound of opening doors and sound of cars. Before we test and processing signals, we will explain the mathematically bases of DWT and HHT transforms and after that we will implement the methods and see the results.In many applications, we need information about time and frequency of a signal. Signal itself does not provide much information and we need to find some methods to determine time-frequency representation of signal.In mathematics, there are methods to describe a function (or signal) in the terms of orthogonal bases functions. Fourier transforms (FT), wavelet transforms (WT), Hilbert-Huang transform (HHT) are examples of those methods to represent and approximate functions. Fourier transforms (FT) and wavelet transforms (WT) are based on Hilbert space and inner product. These tools have been used for many years. The bases function in WT and FT are complete, orthonormal bases or overcomplete.[see Appendix A]Wavelet has many applications and especially the Discrete Wavelet Transform (DWT) is commonly used in engineering, computer science and scientific research. Generally DWT is used for data compression; one example of using DWT is JPEG 2000 which is image compression.The software implementation of DWT is easy, especially in the matlab. One can also implement DWT easily in other programs, such C++ or python.In the many other areas like: seismic geophysics, general signal processes, speech recognition, computer graphic, etc. DWT have been powerful tool to analyse and process data.The Hilbert-Huang transform (HHT) is an empirical method to analyse data. The traditionally methods like WT and FT can only apply to linear processes and has been used in many years for processing data. But question is what can we do about non-linear processes?Almost all processes in the nature are non linear and non stationary, and to analyze this kind of signals, we need other approaches. The need of methods, which can apply to non linear processes, is necessary. HHT is one of those methods, which has been used in recent years. This method has ability to processing both non linear and non stationary processes and that shows the importance and power of HHT.HHT has very simple algorithm to understand and this algorithm can be applied easily to analyze signals. However there are a lot of works to be done. Most work with the HHT was related to its application and there is not so much work about mathematical and theoretical foundation of HHT.Table 1 presents a summary of methods and their comparison:Table 1: Fourier, Wavelet and Hilbert-Huang Transforms Fourier Wavelet Hilbert Basis a priori a priori Adaptive frequency convolution: global uncertainty convolution: regional uncertainty differentiation:local uncertaintyPresentation energy-frequency energy-time-frequency energy-time-frequency Nonlinear no no yes Nonstationary no yes yes Feature Extraction no discrete :no continuous : yesyesTheoretical base theory complete theory complete empiricalFirst the Fourier transform, DWT and HHT will be explained in this chapter and after that in chapter 2, the methods will be implement. Finally in the chapter 3 the results of the methods will be represented1.1 The Fourier transformThe standard Fourier transform is presented in (1.1.1).(Ff )(ω)=12πdte−i ωtf (t )∫ (1.1.1)This can be considered to be only frequency representation of signal and this method can be applied only to stationary signal. The time and frequency can be obtained by cutting off only a well-localized slice (window) of signal and then take its Fourier transform:(T win f )(ω,t )=dsf (s )g (s −t )e −i ωs ∫ (1.1.2)or(T m ,n win)(f )=dsf (s )g (s −nt o )e −im ωo s ∫ (1.1.3)(1.1.2) is the windowed Fourier transform ,which is a standard technique for time frequency localization. Formula (1.1.2) is obtained from (1.1.1) by replacing discrete value (t into nt ).By replacing values: o o ωωm nt t ==, equations’ (1.1.3) can be obtained from (1.1.2). The windowed Fourier transform calls also Short Fast Fourier Transform (SFFT).1.2 Wavelet transformBoth the windowed Fourier transform and wavelet transform give us information in time and frequency domains. The wavelets transform acts similarly to SFFT, but instead of non-local functions, it uses spatially localized functions called wavelets.The continuous wavelet transform formula are shown in (1.2.1) and (1.2.2):(Tf )(a ,b )=a−1/2dtf (t )ψ(t −ba)∫ (1.2.1) Or(T m ,n f )=a o −m /2dtf (t )ψ(a o −m t −nb o )∫(1.2.2) Here b is translation and a is scale. And the function ψ called “mother wavelet” mustsatisfies: dt ψ(t )=0∫ (1.2.3)Here again formula (1.2.2) is obtained from (1.2.1) by replacing discrete value: a =a o m ,b =nb o a o m and n,m belongs Z.The (1.2.3) shows orthogonally of bases functions.We generate b a ,ψ, which is a doubly-index family of wavelet from ψ by scaling (by a) and translating (by b):ψa ,b (x )=a−1/2ψ(x −ba) (1.2.4) a,b ∈R and a ≠0.[1]A function f ∈L 2(ℜ) and its continuous wavelet transform with respect to this wavelet family (1.2.4) are:(T wav f )(a ,b )=f ,ψa ,b =dxf (x )a−1/2ψ(x −ba)∫ (1.2.5)Where ψ is complex conjugate of ψ.The wavelet analyze must satisfy the admissibility condition, namelyC ψ=2πd ζζ−1ˆ ψ (ζ)2<∞∫ (1.2.6)ˆ ψis Fourier transform of ψ. This condition is very important when we want to reconstruct original signal from its transformed .The condition also guarantees theresolution of identity.To recovering function f (original signal), one can use the formula (1.2.7).f =C ψ−1da a 2∫dadb a 2ˆ f (a ,b )a −1/2ψ(t −a b)∫ (1.2.7)Which is inverse transform of ˆ f. Now it becomes clear the importance of condition (1.2.6) i.e. 0<C ψ<∞The definition of parameters a and b is:•Scaling parameter (a=1/f ) describe the length of wavelet, the frequency information of signal.Scaling parameter with high frequency compresses the signal and gives theapproximation of signal and with low frequency, expand signal and providedetailed information.•Translating parameter (b) is related to position and give us information about time in the wavelet transform.Both Fourier transform and wavelet transform take inner product of f with the bases functions g and ψ respectively.There are many mother wavelets and some of them are shown in figure 1.3.1:Meyer MorletHatHaarMexicansFigure 1.3.1 Wavelet basesIn this work Haar wavelet transform is used because it is easy and suitable for signals, and it is easy to implement in hardware.Haar function is defined in equation (1.3.9) and figure 1.3.2 shows the Haar functionψ(x )=10≤x <12−112≤x <10otherwise⎧ ⎨ ⎪ ⎪ ⎩⎪ ⎪ (1.3.9)Figure 1.3.2 Haar functionIn wavelet analyze, any base function (mother wavelet))(,t b a ψ, has to be: • orthonormal,• and any function in the Hilbert space,f ∈L 2(ℜ) can approximate by linear combination of )(,t b a ψ with small error.Haar function satisfies both conditions. For more information see Appendix A1.3 Discrete wavelet transformThere many way to calculate DWT of signal. Here we use a series of filter (high filter and low filter). The DWT of a signal can be computed by passing it through a series filters and produce approximation and details coefficients.In DWT analysis, we often speak about approximation and details. The approximation are the high scale, low frequency part of the signal (low frequency content is the most important part for many signals and provides identity of signal ). The details are low scale, high frequency part (the noise part of signal).The first DWT, was invented by Hungarian mathematician Alfred Haar(1909). The DWT was formulated by the Belgian mathematician Ingrid Daubechies 1988.Daubechies derived a family of wavelets, the first which was the Haar wavelet. There are some wavelet families named Daubechies wavelets after her name. She laid the mathematical foundation to the wavelet transform and explained all mathematically aspects of the DWT.There are two kind of filter being used .The signal can be decomposed by low passing filter g[n] ,and high passing filter h[n]:y low n []=(x ∗g )n []=x k []g n −k []k =−∞∞∑y high n []=(x ∗h )n []=x k []h n −k []k =−∞∞∑The outputs give, the detail coefficients (from the high-pass filter) and approximation coefficients (from low-pass filter). Since half the frequencies of the signal are filteredaway, half of sample can be vanished according to Nyquist’s rule, and then outputs are multiplied by 2:y low n []=(x ∗g )n []=x k []g 2n −k []k =−∞∞∑y high n []=(x ∗h )n []=x k []h 2n −k []k =−∞∞∑1.3.1 One-stage filteringThe original signal passes through two filter and generates two signals.If we actually perform this process on real signal, we get twice as much data as we have in the original signal. Suppose, the original signal has 1000 data, then outputs data, will each have 1000 data. But if we keep only one point out of two in each of the two 2000-data , then we can get the complete information .It produces twosequences called cA (approximation coefficient) and cD (approximation coefficient) in MatLab.1.3.2 Multi-Level DecompositionThe decomposition can be continued and the approximations can bedecomposed in the second level into approximation and detail , and process repeats, so the signal is decomposed into many lower resolution components .This process is called multi-level decomposition.We can choose the number of level based on entropy. The entropy-based criteria determine whether a new decomposition is necessary to obtain minimum-entropy decomposition.The entropy is a common concept used in many fields, mainly in signal processing. There are different entropy criteria and here only one is being used. The entropy E must be additive cost function such that E(0)=0.E(s)=E(si)i∑The concentration in l p norm with 1≤pE2(si)=sipsoE2(s)=sii∑2=spp1.3.3 Discrete Wavelet Packet Analysis (DWPA)The wavelet packet analysis is a generalization of wavelet decomposition, which provides a richer range of possibility for signal analysis.In multi-level decomposition only approximation decomposes at each level. In wavelet packet analysis, the detail as well as the approximation can be decomposed into approximation and detail.The wavelet packet decomposition is shown in figure 1.3.8:Figure1.3.8: Wavelet packetThere are 22n −1different ways to decomposition the signal and we have muchdecomposition to select. To select best tree (best possibility), we used entropy-based criterion. In other words we look at each node of decomposition tree and compute the entropy and this entropy has to be less than previous entropy otherwise stop.[3] The figure 1.3.4 shows the calculation of entropy in nodes according entropy-based criterion.Figure 1.3.4 entropy of tree1.3.4 The Hilbert-Huang TransformTo analyze a signal and find its time-frequency representation by HHT two stages are performed.First we define the IMF. The definition of IMF (which is a dataset)is :1) The number of extreme (maximum and minimum points of dataset) and the number of the zero crossing of function must be equal or differ at most by one.2) The mean value of the envelopes defined by local extreme must be zero.The first stage is to calculate and compute the IMF’s of the signal (The EMD process). The EMD process itself can be done in several steps:1. Find the location of all maxima,)(max t x , and minima )(min t x of )(t x .2. Connect all )(max t x pointes with cubic spline (upper envelope) and do the same procedure to )(min t x to produce lower envelope.3. Compute the mean of spline curves at each point:2)()()(min max t x t x t m +=4. Remove the trend, m (t ). Let )()()(t m t x t d −=.5. Now if d(t) is an IMF according the definition of IMF. Let c i (t )=d (t ) and advance i by 1 .Find the residual )()()(t d t x t r −=.If d (t )is not an IMF (according to definition of IMF), Further sifting is needed6. Reapt steps 1 to 5 ()()(t d t x ←)The sifting process can be stopped when r i becomes monotonic function and no more IMF’s can be extracted.Now we have several IMF (IMFs represented by c i ):nni i nn n r c t x thenc r r c r r c t x r +=−=−=−=∑=−1121211)()(MThe second stage of the HHT process is computing the amplitude and frequency from each IMF’s. And that follows:1. Calculate DFT (discrete Fourier transform) of IMF’s.Then we define Hilbert Transform (HT) and Like others transforms, the Hilbert transform is convolution of function IMF ()()(t c t IMF i =)with g (t )=1t. i.e. []∫∞∞−−==Ητττπd t c PVt g c t c i i i )(1))(*()(2. Calculate the HT. The real and imaginary parts of step1’s DFT consider ascoefficient (M=N/2).y j =1N (Re(X k 1sin(2πk 1j N ))+Im(k 1=0M∑X k 1cos(2πk 1j N))+−1N (Re(X k 2sin(2πk 2j N +Im(k 2=M +1N −1∑X k 2cos(2πk 2j N))3. Calculate the phase φj=tan−1(yjxj) and x j ,y j are from complex numberzj =xj+iyjThe frequency can calculated by taking the derivate of the phasef j =12πdφjdtand amplitude is aj =xj2+yj2 .[1]2. ImplementationIn this part we present the implementation of DWT, DWPA, HHT and NN using MatLab and HHT-program from NASA. In all implementation of DWT Haar wavelet are used. All simulations and computations for DWT are done in MatLab.For HHT (Hilbert-Huang transform), we use HHT program, developed by NASA. Working and analyzing signals with this program is very easy and one can learn it quickly. The HHT is a graphical interface program and one can handle it. The only things you need to do, is loading signals then analyze them. HHT offers all solution to signals processing needs.2.1 Implementation of Wavelet (DWT and DWPA)When DWT and DWPA are applied to the signals, it is obvious that, wavelet packet provides a more flexible tool to the time-frequency analysis of signals. In the general analysis with DWPA is richer than DWT. DWT is actually a subset of the wavelet packet transform (DWPA).There are some disadvantages with DWPA. In DWPA you have many chooses of decomposition and selecting optimal decomposition takes time, compared with DWT. Accuracy and affectivity improves in the DWPA but not so much when you are using DWPA in low levels. We will see those properties by analyzing the real signals. Regarding the choosing number of level of decomposition, we use an entropy-based criterion to select the most optimal decomposition. But there are vary algorithms and criterion to choose best tree. However in the low level decomposition of signals, as we said before using DWPA is not preferable.2.1.1 MethodIn all implementation of DWT and DWPA MatLab is used. We write m files which are shown in figure …(flow diagram) and the codes; can be found in appendix C.We test the codes with noisy sinus signal and after with infrasound signal. From our supervisor we get 6 signals representing two different events see table 2.Table 2. Infrasound signals Signal 1 Signal 2 Signal 3 Event 1 a1 a2 a2 Event 2 b1 b2 b3The program works in the several steps, for example using signal a1 in file a1.m:• Load signal file, here a1.m (all data files must have .m format or you can import data to the workspace).• Chose level of decomposing (for example 3 level) and wavelet base (Haar wavelet).• Decompose the signal and you get details and approximation coefficient. • Threshold the coefficient (kill coefficient which are lower than threshold). Regarding the threshold number se appendix C• Reconstruct the signal which has lower non-zero coefficient than original signal• Compare the energy (norm of signal) of reconstruct signal with original signal and it has to be more than 90% (you can chose any number, but 90% is good enough) of original signal if not go to step 2 and repeat all process (chose a higher level).• Now you have a signal which maintain most energy and has low number non-zero coefficient, in other word you compress the signal without losing so much significant information.Now we have coefficients vector for representing the signal a1.m and we save this vector in as A vector. The same processes are done for others signals and we save their coefficients in new vectors for later applications.In the DWPA, we get quite similar results and we do this in the same program (Matlab program in Appendix C) and you can check result for both methods in the next chapter.2.2 Implementation of Hilbert-Huang transformThe empirical mode data decompositions (EMD) and Hilbert-Huang transform is a general time-frequency method, which consists of two stages in analyzing time series.For a given signal like a1.m, EMD ends up with a representation of the form:y (t )=c i i =1n∑+r nwhere c i is the modes (with zero mean amplitude and is same as IMF i ) and r n is residual.There is a graphical user interface (GUI) program from NASA to analyze signals based on HHT. We use this program to generate all IMFs for each signal.Now we can load signal and decompose signal. Every signal can be decomposed and represented by summation of IMF’s plus residual (r n) and one can collect all IMF’s in the matrix (IMF i rows of matrix and call it for example a1_IMF.m for a1.m signal). Every signal is represented by its matrix, but this matrix has different dimension (like 9x2048 or 7x2048) depending the nature of signal. All signals have its IMF’s matrix and we save them in.It is important to note that matrixes have high dimension (like 9x2048) and it is need to reduce the dimension of matrixes without losing significant information of original signals.Principal component Analysis can be used to reduce the dimension of data set before processing the data with neural networks and it will be explain in the next section. 2.3 Principal Components (PCA)The main idea behind using PCA is to reduce the dimensionality of set of data. This is particularly advantageous if a data set with many variables lies, in reality, close to a two-dimensional sub space (plans).In the case the data can be plotted with respect to these two dimensions, instead of appearing as a large mass numbers (for more information and matlab codes se Appendix B).In many applications, the first two PC’s for the present data account mostly of the total variation, so that a two-dimensional plot with respect to these PC’s, representing most (about 80% , depend to data) of the variation and this gives a reasonably good approximation to the relative positions of the observations in high-dimensional space.In our case we do not need to program some code because PCA is implemented automatically as pre-processing in NN tools in MatLab.2.4 Neural NetworksWe use the interface (GUI) in MatLab … Här skall du skriva enkel hur du använder NN in MatLabsignal a2signal a3signal b1signal b2signal b33. ResultsThe signals of infrasound are shown in figure 3.1 and they are six signals of two different sources.Figure3.1Infrasound signals of eventsWe plot the original signal and reconstructed signal without noise. Notice how it has removed the noise without losing the details and energy (see Fig 3.2, 3.2) WPT gives a better representation of signal but it costs more time (and more operations and more memory to carry out ) see table 1 & 2.Figure 3.2 shows the original signal a1 and its de-noised signal with DWT and PWT. In the figure 3.3 the original signal and its de-noised signal with both methods namely DWT and PWT are shown.Original signal a1Original signal a1Figure 3.2 Decomposition of signal (a1) with discrete wavelet analysis (DWT) and wavelet packet (PWT).Original signal b1compressed signal b1, using wavelet(DWT)Original signal b1compressed signal b1, using wavelet packet(PWT)Figure 3.3 Decomposition of signal (b1) with discrete wavelet analysis (DWT) and wavelet packet (PWT).In the figures 3.4 you can see the tree of decompositions of signal a1 and this tree is optimal tree of decomposition, and in the right side you can se the plot of data in node number 7 or node (3,0), actually this is the plot of coefficient of signal.a1.Figure 3.4 Wavelet packets organized in the tree; Wij (i defines depth and j defines position in the tree) for signal a1.In the figures 3.5 you see optimal tree of decompositions of signal b1 and right side is the plot of data in node number 7 or node (3,0),and these data are coefficient of signal.b1Figure 3.5 Wavelet packets organized in the tree; for signal b1.If WPA and DWT apply to signal b1, gives us the same result, so WPT is not preferable .Actually the DWT, faster and better than WPT for decomposition of signals in the some cases. Choosing which method is optimal is difficult and it depend the nature signals. However the WPT has all advantages of DWT. (se table 1, 2).Table 1 Comparison between DET and WPT at level N=3 for signal a1 Number of zero coefficients percentage L2 recovery (energy) percentage Time (second) Table 2 Comparison between DET and WPT at level N=3 for signal b1 Number of zero coefficients percentage L2 recovery (energy) percentage Time (second)DWT 89.7949 97.0663 0.0801WPT 91.6992 97.0902 0.159536DWT 75.9766 99.1707 0.031624WPT 95.9766 99.1707 0.148443Hilbert-Huang transformationThe results of decomposing signals into theirs IMF’s plots , with HHT are shown in Figures 2.2.1 and 2.2.2 is plot of IMF’s.0.1 0 -0.1 0.05 0 -0.05 0.05 0 -0.05 0.05 0 -0.05 0.2 0 -0.2 0.2 0 -0.2 0.05 0 -0.05 0.05 0 -0.05 0.01 0 -0.01 0.02 0 -0.02 1 0 -1 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500Figure 2.2.1:IMF’s of the signal for signal a11 0 -1 1 0 -1 0.5 0 -0.5 0.2 0 -0.2 0.2 0 -0.2 0.1 0 -0.1 0.02 0 -0.02 -0.02 -0.03 -0.04 1 0 -1 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500Figure 2.2.2:IMF’s of the signal for signal b1Graphical Representation of Data using Principal ComponentsA plot of these data with respect to the first PC for sex signals was given in Figure 2.3.1siganl a1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -0.5 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.5siganl a2 1siganl a30.5 Second Principal Component 0 0.5 1 1.5 First Principal Component 2 2.5Second Principal ComponentSecond Principal Component0-0.5-1-1.5-200.5 1 1.5 First Principal Component2-2.5 -101 2 3 First Principal Component4siganl b1 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -1 0.8 0.6 0.4 Second Principal Component 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 1 2 3 First Principal Component 4 -1.2 -2siganl b2 50 40 30 Second Principal Component 20 10 0 -10 -20 -30 -40 -1 0 1 First Principal Component 2 -50 -50siganl b3Second Principal Component0 50 First Principal Component100For more information about PCA and matlab code se Appendix B & DConclusionThe HHT method is very powerful tool to processing signals and the HHT Algorithm can easily be implemented in any programming languages like C++ or matlab itself. However the mathematical foundation of HHT is not complete but we can still use it. DWT also provide a time-frequency representations of signals. Using DWT especially in MatLab is practically and easy. Both DWT(discrete wavelet transform) or DWPA (Wavelet Packet Analysis) are using the WT methods and in the some cases we got the same results (even DWT gives better result). There are some similarity between HHT and DWPA in the terms of filter bank and level of decomposition.Appendix A Review of Hilbert space analysisTo understand the wavelet transform, it is helpful to have basic knowledge about Hilbert space analysis. 1. The inner product space. Let V be the space of real-value functions with domain the real line (ℜ) and let the inner product f , g be defined for all f and g in V as:f ,g =∫f ( x ) g( x)dxℜWhere g( x ) is complex conjugate of g(x). f and g are orthogonal if f ,g =0;. 2. Hilbert space of interest 1/ 2 The norm of function f defined as f ≡ f , f and this is analogue to a length in Euclidean space. p Using the norm, a Hilbert space L (V ) can be defined as:L p (V ) ≡ {f ∈ V | fp< ∞}If p=2 then the space of interest contains all f that are square-integrable on 2 V. We concentrate on the Hilbert space L (ℜ) . 3. Specify a linear operator on the Hilbert space A linear operator L in Hilbert space H is one for whichL( f + g) = L( f ) + L(g) and L(cf ) = cL( f )for all constants c and all f,g in H. The eigenvalues of a linear operator L are generally found by solving the equationLf = λfy ∈ Dom(L)to determine the allowed Eigen values λ1 , λ 2 ,..., λ n and then finding the associated eigenfunctions ϕ1 ,ϕ 2 ,...,ϕ n which satisfy the above equation.Wavelet analysis techniquesIn Hilbert space analysis, one finds the union of eigenspaces of a linear the space operator and determines the basis of that space. The operator itself is chosen to。
同步压缩小波变换matlab程序英文回答:Wavelet transform is a powerful tool in signal processing and data compression. It is widely used in various fields such as image and audio compression, denoising, and feature extraction. In MATLAB, there are built-in functions and toolboxes that can be used to perform wavelet transform and compression.To perform synchronous wavelet compression in MATLAB, we can follow these steps:1. Load the signal or image data: We first need to load the signal or image data that we want to compress. This can be done using the appropriate MATLAB functions, such as`audioread` for audio signals or `imread` for images.2. Choose a wavelet: Next, we need to choose a suitable wavelet for the compression. MATLAB provides a variety ofwavelets, such as Daubechies, Coiflets, and Symlets. We can use the `wfilters` function to obtain the coefficients of a specific wavelet.3. Perform wavelet decomposition: We can use the`wavedec` function to decompose the signal or image into different frequency subbands using the chosen wavelet. This will result in a set of approximation and detail coefficients.4. Set a compression threshold: In order to reduce the amount of data to be stored or transmitted, we can set a compression threshold to discard or truncate the detail coefficients with small magnitudes. This can be done by comparing the magnitude of each coefficient with the threshold value.5. Reconstruct the compressed signal or image: After discarding or truncating the detail coefficients, we can use the `waverec` function to reconstruct the compressed signal or image using the remaining approximation anddetail coefficients.6. Evaluate the compression performance: Finally, wecan evaluate the compression performance by comparing the quality of the reconstructed signal or image with the original data. This can be done using various metrics such as peak signal-to-noise ratio (PSNR) or mean squared error (MSE).中文回答:小波变换是信号处理和数据压缩中的一种强大工具。