小波包分解
小波分析的语音信号噪声消除方法
小波分析的语音信号噪声消除方法小波分析是一种有效的信号处理方法，可以用于噪声消除。
在语音信号处理中，噪声常常会影响语音信号的质量和可理解性，因此消除噪声对于语音信号的处理非常重要。
下面将介绍几种利用小波分析的语音信号噪声消除方法。
一、阈值方法阈值方法是一种简单而有效的噪声消除方法，它基于小波变换将语音信号分解为多个频带，然后通过设置阈值将各个频带的噪声成分消除。
1.1离散小波变换（DWT）首先，对语音信号进行离散小波变换（DWT），将信号分解为近似系数和细节系数。
近似系数包含信号的低频成分，而细节系数包含信号的高频成分和噪声。
1.2设置阈值对细节系数进行阈值处理，将细节系数中幅值低于设定阈值的部分置零。
这样可以将噪声成分消除，同时保留声音信号的特征。
1.3逆变换将处理后的系数进行逆变换，得到去噪后的语音信号。
1.4优化阈值选择为了提高去噪效果，可以通过优化阈值选择方法来确定最佳的阈值。
常见的选择方法有软阈值和硬阈值。
1.4.1软阈值软阈值将细节系数进行映射，对于小于阈值的细节系数，将其幅值缩小到零。
这样可以在抑制噪声的同时保留语音信号的细节。
1.4.2硬阈值硬阈值将细节系数进行二值化处理，对于小于阈值的细节系数，将其置零。
这样可以更彻底地消除噪声，但可能会损失一些语音信号的细节。
二、小波包变换小波包变换是对离散小波变换的改进和扩展，可以提供更好的频带分析。
在语音信号噪声消除中，小波包变换可以用于更精细的频带选择和噪声消除。
2.1小波包分解将语音信号进行小波包分解，得到多层的近似系数和细节系数。
2.2频带选择根据噪声和语音信号在不同频带上的能量分布特性，选择合适的频带对语音信号进行噪声消除。
2.3阈值处理对选定的频带进行阈值处理，将噪声成分消除。
2.4逆变换对处理后的系数进行逆变换，得到去噪后的语音信号。
三、小波域滤波小波域滤波是一种基于小波变换的滤波方法，通过选择合适的小波函数和滤波器来实现噪声消除。
Matlab中的小波变换与小波包分析方法详解
Matlab中的小波变换与小波包分析方法详解引言近年来，小波变换在信号处理领域中得到了广泛的应用。
小波变换是一种能够捕捉信号时频特性的有效工具，可以用来分析、压缩和去噪各种类型的信号。
本文将详细介绍Matlab中的小波变换和小波包分析方法，以帮助读者更好地理解和应用这一强大的信号处理技术。
一、小波变换（Wavelet Transform）小波变换是一种将信号分解成不同尺度的基函数的技术。
与传统的傅里叶变换相比，小波变换具有更好的时频局部化特性。
Matlab中提供了丰富的小波分析工具箱，可以方便地进行小波变换的计算。
1.1 小波基函数小波基函数是小波变换的基础。
不同类型的小波基函数适用于不同类型的信号。
在Matlab中，我们可以使用多种小波基函数，如Daubechies小波、Haar小波和Morlet小波等。
1.2 小波分解小波分解是指将信号分解成多个尺度的小波系数。
通过小波分解，我们可以获取信号在不同尺度上的时频特性。
Matlab中提供了方便的小波分解函数，例如'dwt'和'wavedec'。
1.3 小波重构小波重构是指根据小波系数重新构建原始信号。
通过小波重构，我们可以恢复原始信号的时域特性。
在Matlab中，可以使用'idwt'和'waverec'函数进行小波重构。
二、小波包分析（Wavelet Packet Analysis）小波包分析是对小波变换的进一步扩展，它允许对信号进行更精细的分解和重构。
小波包分析提供了一种更灵活的信号分析方法，能够获得更详细的时频特性。
2.1 小波包分解小波包分解是指将信号分解成具有不同频带的小波包系数。
与小波分解相比，小波包分解提供了更高的分辨率和更详细的频谱信息。
在Matlab中，可以使用'wavedec'函数进行小波包分解。
2.2 小波包重构小波包重构是根据小波包系数重新构建原始信号。
小波包PPT课件
引言
小波分解示意图----每层分解只对低频部分细分
S
A1
D1
A2
D2
A3
D3
4
引言
小波包分解，在小波分解的基础上进一步细分高频部分，达 到更优的时频局部化效果
S
A1
D1
A2,1
D2,1
A2,2
D2,2
5
A3,1
D3,1
A3,2
D3,2
A3,3
D3,3
A3,4
D3,4
小波包原理
❖ 所谓小波包，简单地说就是一个函数族。由 它们构造出的规范正交基库。从此库中可以 选出的许多规范正交基，小波正交基只是其 中的一组，所以小波包是小波概念的推广。
包，称为小波包系数。G，H为小波分解滤波器， H与尺度函数 有关，G与 j (t)有关。二进小波包 分解的快速算法为：
p01 (t) p 2i 1
j
f
(t) H (k
2t
)
p
i j
1
(t
)
k
p
2i j
k
G(k
2t
)
p
i j
1
(t
)
9
重构算法为：
p
i j
(t
)
2[
h(t
2k
)
p
2 i 1 j 1
(t
)
g
(t
2k
)
p
2i j 1
(t
)]
k
k
式中，j J 1, J 2,,1，0；i 2 j ,2 j1,,2,1;
J
log
N 2
, h,
g为小波重构滤波器，
基于小波包分解的整体变分去噪算法
Ｃ ｏ ｍｐ ｕ ｔ ｅ ｒ Ｅｎ ｇ ｉ ｎ ｅ ｅ ｒ ｉ ｎ ｇ ａ ｎ ｄ Ａｐ ｐ ｌ ｉ ｃ ａ ｔ ｉ ｏ ｎ ｓ ， ２ ０ １ ３ ， ４ ９ （ ６ ） ： １ ５ ６ — １ ５ ８ ．
Ａ ｂ ｓ ｔ ｒ ａ ｃ ｔ ： Ｔ ｈ ｅ Ｔ ｏ ｔ ａ ｌ Ｖａ ｒ ｉ ￣ｉ ｏ ｎ （ Ｔ Ｖ） ｍｏ ｄ ｅ ｌ ｏ ｆ Ｒ ｕ ｄ ｉ ｎ ｅ ｔ ａ １ ． ｆ ｏ ｒ ｉ ｍａ ｇ ｅ ｄ ｅ ｎ ｏ ｉ ｓ ｉ ｎ ｇ ｉ ｓ ｃ ｏ ｎ ｓ ｉ ｄ ｅ ｒ ｅ ｄ ａ ｓ ｏ ｎ ｅ ｏ ｆ ｔ ｈ ｅ ｂ ｅ ｓ ｔ ｄ ｅ ｎ ｏ ｉ ｓ ｉ ｎ ｇ ｍｏ ｄ ｅ ｌ ｓ ．
西 北工 业 大学 理学 院 数学 系 ， 西安 ７ １ ０ １ ２ ９
Ｄｅ ｐ ａ ｒ ｔ ｍｅ ｎ ｔ ｏ ｆ Ａｐ ｐ ｌ ｉ ｅ ｄ Ｍａ ｔ ｈ ｅ ｍａ ｔ ｉ ｃ ｓ ， Ｃｏ ｌ ｌ ｅ ｇ ｅ ｏ ｆ Ｓ ｃ ｉ ｅ ｎ ｃ ｅ ， Ｎｏ ｒ ｔ ｈ ｗｅ ｓ ｔ ｅ ｒ ｎ Ｐ ｏ ｌ ｙ ｔ ｅ ｃ ｈ ｎ ｉ ｃ ａ ｌ Ｕｎ ｉ ｖ ｅ ｒ ｓ ｉ ｔ ｙ ， Ｘｉ ’ ａ ｎ ７ １ ０１ ２ ９ ， Ｃｈ ｉ ｎ ａ
Ｔ ＡＮＧ Ｃｈａ ｎ ｇ ｌ ｉ ｎｇ ，Ｐ ＥＮＧ Ｇｕ ｏ ｈ ｕａ ．Ｔ ｏ ｔ ａ ｌ ｖ ａ ｒ ｉ ａ ｔ ｉ ｏ ｎ ａ ｌ ｇ ｏ ｒ ｉ ｔ ｈ ｍ ｏ ｆ ｉ ｍａ ｇ ｅ ｄ ｅ ｎ ｏ ｉ ｓ ｉ ｎ ｇ ｂ ａ ｓ ｅ ｄ ｏ ｎ ｗａ ｖ ｅ ｌ ｅ ｔ ｐａ ｃ ｋｅ ｔ ｄ ｅ ｃ ｏ ｍｐ ｏ ｓ ｉ ｔ ｉ ｏ ｎ．
小波分析及小波包分析
小波分析及小波包分析在利用matlab做小波分析时，小波分解函数和系数提取函数的结果都是分解系数。
我们知道，复杂的周期信号可以分解为一组正弦函数之和，及傅里叶级数，而傅里叶变换对应于傅里叶级数的系数；同样，信号也可以表示为一组小波基函数之和，小波变换系数对应于这组小波基函数的系数。
多尺度分解是按照多分辨分析理论，分解尺度越大，分解系数的长度越小（是上一个尺度的二分之一）。
我们会发现分解得到的小波低频系数的变化规律和原始信号相似，但要注意低频系数的数值和长度与原始信号以及后面重构得到的各层信号是不一样的。
小波分解：具体实现过程可以分别设计高通滤波器和低通滤波器，得到高频系数和低频系数，并且每分解一次数据的长度减半。
小波重构，为分分解的逆过程，先进行增采样，及在每两个数之间插入一个0，与共轭滤波器卷积，最后对卷积结果求和。
在应用程中，我们经常利用各层系数对信号进行重构（注意虽然系数数少于原信号点数，但是重构后的长度是一样的），从而可以有选择的观看每一频段的时域波形。
从而确定冲击成分所在频率范围。
便于更直观的理解，小波分解，利用各层系数进行信号重构过程我们可以认为是将信号通过一系列的不同类型的滤波器，从而得到不同频率范围内的信号，及将信号分解。
小波消噪:运用小波分析进行一维信号消噪处理和压缩处理，是小波分析的两个重要的应用。
使用小波分析可以将原始信号分解为一系列的近似分量和细节分量，信号的噪声主要集中表现在信号的细节分量上。
使用一定的阈值处理细节分量后，再经过小波重构就可以得到平滑的信号。
小波常用函数[C,L]=wavedec(s,3,'db1');%用小波函数db1对信号s进行3尺度分解其中C为分解后低频和高频系数，L存储低频和高频系数的长度。
X=wrcoef(‘type’,C,L,’wname’,N)%对一维小波系数进行单支重构，其中N表示对第几层的小波进行重构X=wrcoef(‘a’,C,L,’wname’,3)%对第三层的低频信号进行重构，如果a变为d的话，表示对低频分量进行重构。
小波包变换和小波变换
小波包变换和小波变换小波包变换和小波变换是一种信号分析和处理的方法，它们可以将信号分解成不同尺度和频率的成分，并可以分析和处理这些成分。
下面将对小波包变换和小波变换进行解释。
1. 小波包变换：小波包变换是在小波变换的基础上发展而来的一种方法。
小波包变换将信号分解成多个子带，并对每个子带进行进一步的分解。
相比于小波变换，小波包变换提供了更高的频率分辨率和更细的频率划分。
小波包变换的核心思想是使用不同的小波基函数对信号进行分解。
通过选择不同的小波基函数，可以获得不同尺度和频率的信号成分。
小波包变换可以通过反复迭代的方式，不断将信号分解成更细的频率带，进一步提高频率分辨率。
在每一级分解中，信号被分解成低频和高频两部分，低频部分可以继续进行进一步的分解。
小波包变换的优势在于能够提供更详细的频域信息，可以更好地分析信号的特征和结构。
它在信号处理、图像处理等领域有着广泛的应用，例如信号去噪、特征提取等。
2. 小波变换：小波变换是一种将信号分解成不同频率成分的方法。
通过小波变换，我们可以将信号从时域转换到频域，同时可以分析信号的时间和频率特性。
小波变换的基本思想是使用小波基函数对信号进行分解。
小波基函数是一种具有局部性质的函数，它能够在时域和频域中同时提供较好的分辨率。
通过选择不同的小波基函数，可以获得不同频率范围内的信号成分。
小波变换通过对信号进行连续的分解和重构，可以分析信号的频域特性。
小波变换有多种变体，其中最常用的是离散小波变换（DWT）。
离散小波变换将信号分解成多个尺度和频率的子带，通过这些子带可以分析信号的不同频率成分。
离散小波变换具有高效性和局部性，可以在信号处理中广泛应用，例如信号去噪、压缩等。
总结：小波包变换是在小波变换的基础上发展的一种方法，它能够提供更高的频率分辨率和更细的频率划分。
小波包变换通过选择不同的小波基函数，将信号分解成多个子带，并对每个子带进行进一步的分解。
相比之下，小波变换是将信号分解成不同频率成分的方法，通过选择不同的小波基函数，可以获得不同频率范围内的信号成分。
小波包变换的基本原理和使用方法
小波包变换的基本原理和使用方法引言：小波包变换（Wavelet Packet Transform）是一种信号分析技术，它在小波变换的基础上进一步拓展，能够提供更丰富的频域和时域信息。
本文将介绍小波包变换的基本原理和使用方法，帮助读者更好地理解和应用这一技术。
一、小波包变换的基本原理小波包变换是一种多分辨率分析方法，它利用小波基函数对信号进行分解和重构。
与传统的傅里叶变换相比，小波包变换能够提供更精细的频域和时域信息，适用于非平稳信号的分析。
小波包变换的基本原理如下：1. 信号分解：首先将原始信号分解为不同频率的子信号，通过迭代地将信号分解为低频和高频部分，形成小波包树结构。
2. 小波基函数：在每一层分解中，选取合适的小波基函数进行信号分解。
小波基函数具有局部性和多分辨率特性，能够更好地捕捉信号的局部特征。
3. 分解系数：分解过程中，每个子信号都会生成一组分解系数，用于表示信号在不同频率上的能量分布。
分解系数可以通过滤波和下采样得到。
二、小波包变换的使用方法小波包变换在信号处理领域有广泛的应用，包括信号去噪、特征提取、模式识别等。
下面将介绍小波包变换的常见使用方法。
1. 信号去噪：小波包变换可以提供更丰富的频域和时域信息，因此在信号去噪领域有较好的效果。
通过对信号进行小波包分解，可以将噪声和信号分离，然后对噪声进行滤波处理，最后通过重构得到去噪后的信号。
2. 特征提取：小波包变换可以提取信号的局部特征，对于信号的频率变化和时域特征有较好的描述能力。
通过对信号进行小波包分解，可以得到不同频率下的分解系数，进而提取出信号的主要特征。
3. 模式识别：小波包变换在模式识别中也有广泛的应用。
通过对信号进行小波包分解，可以得到不同频率下的分解系数，进而提取出信号的特征向量。
利用这些特征向量，可以进行模式分类和识别。
4. 压缩编码：小波包变换可以将信号进行有效的压缩编码。
通过对信号进行小波包分解，可以将信号的主要信息集中在少量的分解系数中，从而实现信号的压缩。
小波包变换（WaveletPacketTransform）的学习笔记
⼩波包变换（WaveletPacketTransform）的学习笔记对于⼀个连续的周期信号，可以将其分解为⼀组频率不同的三⾓函数信号的线性组合，这就是傅⾥叶级数的本质，将信号从时域投影到频域中的不同频段上来完成分解。
当这个周期信号的周期趋近于⽆穷⼤时，傅⾥叶级数就变成了傅⾥叶变换。
此时的信号本质上是⼀个连续⾮周期信号，傅⾥叶变换的意义就在于对其进⾏分解，同样也是以⼀组三⾓函数作为正交基，并通过这组三⾓函数基的线性组合来表⽰原信号。
数学表达为：由于三⾓函数是⼀个⽆限长的信号，在时域上不具有局部性，因此以其作为正交基对信号进⾏拟合时，具有以下两个不⾜：第⼀，对于突变信号，如阶跃信号或尖峰信号，其需要⼤量的三⾓函数基进⾏组合才能完成较好的信号拟合；第⼆，由于三⾓函数不具备在时域上的局部性，因此在对信号进⾏傅⾥叶变换时，仅仅只能获取到信号在频域上的分布信息，并不能获取到这些不同频率的信号分量在时域上出现的位置。
因此傅⾥叶变换对于⾮平稳信号的分解会遗失其在时域上的变化信息。
⼩波变换就是为了解决对⾮平稳信号的分解问题⽽产⽣的数学⽅法。
相⽐于傅⾥叶变换使⽤⼀组⽆限长的三⾓函数基进⾏信号拟合，⼩波变换使⽤的是⼀组正交的、迅速衰减的⼩波函数基进⾏信号拟合。
这种⼩波函数基可通过其尺度变量和平移变量，获得不同的频率和时间位置。
因此在利⽤这种⼩波函数基对信号进⾏分解时，可以⽤较少的⼩波函数基就拟合出突变信号（稀疏编码特性），同时也能获得不同频率的信号分量在时域上的出现位置。
⽤于⽣成⼀组不同频率和时移的⼩波函数的⼩波函数，称为基本⼩波（Basic Wavelet），由其⽣成的⼀组⼩波函数，是该基本⼩波的⼀个⼩波族（Wavelet Family），表⽰为：，其中为尺度参数，通过伸缩控制⼩波的尺度（频率），为平移参数，通过移位控制⼩波在时域中的出现位置。
这两个参数的作⽤顺序是先作平移，再作伸缩。
对这⼀族⼩波函数进⾏归⼀化，即得到⼀组⼩波函数基。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%小波包分解程序%% m=load('300_30.txt'); 导入文件名为300_30的txt文件 N=length(m); for i=1:N-1 ; q(i,1)=m(i,1); end; d=q'; s1=d; change=1000; [c,l] = wavedec(d,3,'db4'); %提取小波分解后的低频系数 ca3=appcoef(c,l,'db4',3); %提取各层小波分解后的高频系数 cd3=detcoef(c,l,3); cd2=detcoef(c,l,2); cd1=detcoef(c,l,1); %对信号强制消噪 cdd3=zeros(1,length(cd3));%第三层高频系数cd3全置0 cdd2=zeros(1,length(cd2));%第二层高频系数cd2全置0 cdd1=zeros(1,length(cd1));%第一层高频系数cd1全置0
c1=[ca3,cdd3,cdd2,cdd1];%建立新的系数矩阵 s2=waverec(c1,l,'db4')%为新的分解结构
%[thr,sorh,keepapp]=ddencmp('den','wv',d); %s2=wdencmp('gbl',c,l,'db4',4,thr,sorh,keepapp); %subplot(413) %plot(1:change,s2(1:change)); %title('默认软阈值消噪后信号')
figure(1) subplot(9,2,1) plot(1:change,s1(1:change)) title('原始信号') ylabel('S1')
subplot(9,2,2) plot(1:change,s2(1:change)) title('强制消噪后信号') ylabel('S2')
wpt=wpdec(s1,3,'db1','shannon'); %plot(wpt); %重构第三层个节点小波系数 s130=wprcoef(wpt,[3,0]); s131=wprcoef(wpt,[3,1]); s132=wprcoef(wpt,[3,2]); s133=wprcoef(wpt,[3,3]); s134=wprcoef(wpt,[3,4]); s135=wprcoef(wpt,[3,5]); s136=wprcoef(wpt,[3,6]); s137=wprcoef(wpt,[3,7]); %计算第三层个节点小波能量 s10=norm(s130); s11=norm(s131); s12=norm(s132); s13=norm(s133); s14=norm(s134); s15=norm(s135); s16=norm(s136); s17=norm(s137); %计算方差 st10=std(s130); st11=std(s131); st12=std(s132); st13=std(s133); st14=std(s134); st15=std(s135); st16=std(s136); st17=std(s137);
disp('正常信号的特征向量'); snorm1=[s10,s11,s12,s13,s14,s15,s16,s17]; std1=[st10,st11,st12,st13,st14,st15,st16,st17]; %显示三层个节点小波系数 subplot(9,2,3);plot(1:change,s130(1:change)); ylabel('S130'); subplot(9,2,5);plot(1:change,s131(1:change)); ylabel('S131'); subplot(9,2,7);plot(1:change,s132(1:change)); ylabel('S132'); subplot(9,2,9);plot(1:change,s133(1:change)); ylabel('S133'); subplot(9,2,11);plot(1:change,s134(1:change)); ylabel('S134'); subplot(9,2,13);plot(1:change,s135(1:change)); ylabel('S135'); subplot(9,2,15);plot(1:change,s136(1:change)); ylabel('S136'); subplot(9,2,17);plot(1:change,s137(1:change)); ylabel('S137');
wpt=wpdec(s2,3,'db1','shannon'); %plot(wpt); s230=wprcoef(wpt,[3,0]); s231=wprcoef(wpt,[3,1]); s232=wprcoef(wpt,[3,2]); s233=wprcoef(wpt,[3,3]); s234=wprcoef(wpt,[3,4]); s235=wprcoef(wpt,[3,5]); s236=wprcoef(wpt,[3,6]); s237=wprcoef(wpt,[3,7]);
s20=norm(s230); s21=norm(s231); s22=norm(s232); s23=norm(s233); s24=norm(s234); s25=norm(s235); s26=norm(s236); s27=norm(s237);
st20=std(s230); st21=std(s231); st22=std(s232); st23=std(s233); st24=std(s234); st25=std(s235); st26=std(s236); st27=std(s237);
disp('故障信号的特征向量'); snorm2=[s20,s21,s22,s23,s24,s25,s26,s27]; std2=[st20,st21,st22,st23,st24,st25,st26,st27];
subplot(9,2,4);plot(1:change,s230(1:change)); ylabel('S230'); subplot(9,2,6);plot(1:change,s231(1:change)); ylabel('S231'); subplot(9,2,8);plot(1:change,s232(1:change)); ylabel('S232'); subplot(9,2,10);plot(1:change,s233(1:change)); ylabel('S233'); subplot(9,2,12);plot(1:change,s234(1:change)); ylabel('S234'); subplot(9,2,14);plot(1:change,s235(1:change)); ylabel('S235'); subplot(9,2,16);plot(1:change,s236(1:change)); ylabel('S236'); subplot(9,2,18);plot(1:change,s237(1:change)); ylabel('S237');
%fft N=1024; figure(2) y1=fft(s1,N); py1=abs(y1) %py1=y1.*conj(y1)/N; y2=fft(s2,N); py2=abs(y2) %py2=y2.*conj(y2)/N;
y130=fft(s130,N); py130=y130.*conj(y130)/N; y131=fft(s131,N); py131=y131.*conj(y131)/N; y132=fft(s132,N); py132=y132.*conj(y132)/N; y133=fft(s133,N); py133=y133.*conj(y133)/N; y134=fft(s134,N); py134=y134.*conj(y134)/N; y135=fft(s135,N); py135=y135.*conj(y135)/N; y136=fft(s136,N); py136=y136.*conj(y136)/N; y137=fft(s137,N); py137=y137.*conj(y137)/N;
y230=fft(s230,N); py230=y230.*conj(y230)/N; y231=fft(s231,N); py231=y231.*conj(y231)/N; y232=fft(s232,N); py232=y232.*conj(y232)/N; y233=fft(s233,N); py233=y233.*conj(y233)/N; y234=fft(s234,N); py234=y234.*conj(y234)/N; y235=fft(s235,N); py235=y235.*conj(y235)/N; y236=fft(s236,N); py236=y236.*conj(y236)/N; y237=fft(s237,N); py237=y237.*conj(y237)/N; para=512; f=1000*(0:para-1))/N; subplot(1,2,1); plot(f,py1(1:para)); ylabel('P1'); title('原始信号的功率谱') subplot(1,2,2); plot(f,py2(1:para)); ylabel('P2'); title('故障信号的功率谱')
figure subplot(4,2,1); plot(f,py130(1:para)); ylabel('P130'); title('S130的功率谱') subplot(4,2,2); plot(f,py131(1:para)); ylabel('P131'); title('S131的功率谱') subplot(4,2,3); plot(f,py132(1:para)); ylabel('P132'); subplot(4,2,4); plot(f,py133(1:para)); ylabel('P133'); subplot(4,2,5); plot(f,py134(1:para)); ylabel('P134'); subplot(4,2,6); plot(f,py135(1:para)); ylabel('P135'); subplot(4,2,7);