急求【MATLAB源代码】 基于小波分析的电力系统谐波检测

t=0:0.0001:1.3

figure

subplot(321)

f1=sin(1*pi*10*t)

plot(f1)

title('频率为5hz的正弦波')

Ylabel('幅值')

Xlabel('时间')

subplot(322)

f2=sin(2*pi*10*t)

plot(f2)

title('频率为10hz的正弦波')

Ylabel('幅值')

Xlabel('时间')

subplot(323)

f3=sin(3*pi*10*t)

plot(f3)

title('频率为15hz的正弦波')

Ylabel('幅值')

Xlabel('时升族弯间')

f=f1+f2+f3

subplot(324)

plot(f)

title('合成的正弦波')

Ylabel('幅值')

Xlabel('时间')

subplot(325)

coefs=cwt(f,[1:1:10],'吵闷db3','plot')

title('对于不同尺度下的小波系数值')

Ylabel('尺穗并度')

Xlabel('时间')

程序是下面这样，但只能处理长宽一样的方形图像，灰度和彩祥宽谨色图像都可，要用其他图像只需把Lena.bmp改为其他图像，但图像要保存在m文件所在路径下。下面巧耐还有一个运行后的图像（之一），希望能对你有所帮助。

clearclc

%%%%%%%%%%测试图像只能是方形图像，长宽像素一样。

f=imread('Lena.bmp')%%读取图像数据，图像只能保存在m文件所在的路径下

d=size(f)

if length(d)>2

f=rgb2gray((f))%%%%%%%%如果是彩色图像则转化为灰度图

end

T=d(1)

SUB_T=T/2

% 2.进行二维小波分解

l=wfilters('db10','l') % db10（消失矩为10)低通分解滤波器冲击响应（长度为20）

L=T-length(l)

l_zeros=[l,zeros(1,L)] % 矩阵行数与输入图像一致，为2的整数幂

h=wfilters('db10','h') % db10（消失矩为10)高通分解滤波器冲击响应（长度为20）

h_zeros=[h,zeros(1,L)] % 矩阵行数与输入图像一致，为2的整数幂

for i=1:T % 列变换

row(1:SUB_T,i)=dyaddown( ifft( fft(l_zeros).*fft(f(:,i)') ) ).' % 圆周卷积<->FFT

row(SUB_T+1:T,i)=dyaddown( ifft( fft(h_zeros).*fft(f(:,i)') ) ).' % 圆周卷积<->FFT

end

for j=1:T % 行变换

line(j,1:SUB_T)=dyaddown( ifft( fft(l_zeros).*fft(row(j,:)) ) ) % 圆周卷积<->FFT

line(j,SUB_T+1:T)=dyaddown( ifft( fft(h_zeros).*fft(row(j,:)) ) ) % 圆周谨基卷积<->FFT

end

decompose_pic=line % 分解矩阵

% 图像分为四块

lt_pic=decompose_pic(1:SUB_T,1:SUB_T) % 在矩阵左上方为低频分量--fi(x)*fi(y)

rt_pic=decompose_pic(1:SUB_T,SUB_T+1:T) % 矩阵右上为--fi(x)*psi(y)

lb_pic=decompose_pic(SUB_T+1:T,1:SUB_T) % 矩阵左下为--psi(x)*fi(y)

rb_pic=decompose_pic(SUB_T+1:T,SUB_T+1:T) % 右下方为高频分量--psi(x)*psi(y)

% 3.分解结果显示

figure(1)

subplot(2,1,1)

imshow(f,[]) % 原始图像

title('original pic')

subplot(2,1,2)

image(abs(decompose_pic)) % 分解后图像

title('decomposed pic')

figure(2)

% colormap(map)

subplot(2,2,1)

imshow(abs(lt_pic),[]) % 左上方为低频分量--fi(x)*fi(y)

title('\Phi(x)*\Phi(y)')

subplot(2,2,2)

imshow(abs(rt_pic),[]) % 矩阵右上为--fi(x)*psi(y)

title('\Phi(x)*\Psi(y)')

subplot(2,2,3)

imshow(abs(lb_pic),[]) % 矩阵左下为--psi(x)*fi(y)

title('\Psi(x)*\Phi(y)')

subplot(2,2,4)

imshow(abs(rb_pic),[]) % 右下方为高频分量--psi(x)*psi(y)

title('\Psi(x)*\Psi(y)')

% 5.重构源图像及结果显示

% construct_pic=decompose_matrix'*decompose_pic*decompose_matrix

l_re=l_zeros(end:-1:1) % 重构低通滤波

l_r=circshift(l_re',1)' % 位置调整

h_re=h_zeros(end:-1:1) % 重构高通滤波

h_r=circshift(h_re',1)' % 位置调整

top_pic=[lt_pic,rt_pic] % 图像上半部分

t=0

for i=1:T % 行插值低频

if (mod(i,2)==0)

topll(i,:)=top_pic(t,:)% 偶数行保持

else

t=t+1

topll(i,:)=zeros(1,T) % 奇数行为零

end

end

for i=1:T % 列变换

topcl_re(:,i)=ifft( fft(l_r).*fft(topll(:,i)') )' % 圆周卷积<->FFT

end

bottom_pic=[lb_pic,rb_pic] % 图像下半部分

t=0

for i=1:T % 行插值高频

if (mod(i,2)==0)

bottomlh(i,:)=bottom_pic(t,:) % 偶数行保持

else

bottomlh(i,:)=zeros(1,T) % 奇数行为零

t=t+1

end

end

这个只是一级分解，matlab自带的函数可以实现多级分解，级数由编程者自己确定。

是的，是一样的。

%含噪声的三角波与正弦波的组合

%利用db5小波对信号进行7层分解

%生产正弦信号

clcclose allclear all

N=1000

t=1:N

sig1=sin(0.3*t)

%生成三角形波形

sig2(1:500)=((1:500)-1)/500

sig2(501:N)=(1000-(501:1000))/500

figure(1)

subplot(211)

plot(t,sig1,'linewidth',2)

xlabel('样本序号 N')

ylabel('幅值A'明灶)

subplot(212)

plot(t,sig2,'linewidth',2)

xlabel('样本序号 N')

ylabel('幅值A')

%叠加信号

x=sig1+sig2+randn(1,N)

figure(2)

plot(t,x,'linewidth',2)

xlabel('样本序号 N')

ylabel('幅含喊值A')%一维小波分解

[c,l]=wavedec(x,7,'db5')%重构第1-7层逼近系数

a7=wrcoef('a',c,l,'db5',7)

a6=wrcoef('a',c,l,'db5',6)

a5=wrcoef('a',c,l,'db5',5)

a4=wrcoef('a',c,l,'db5',4)

a3=wrcoef('a',c,l,'db5',3)

a2=wrcoef('a',c,l,'db5',2)

a1=wrcoef('a',c,l,'db5',1)%显示激老扮逼近系数

figure(3)

subplot(711)

plot(a7,'linewidth',2)

ylabel('a7')

subplot(712)

plot(a6,'linewidth',2)

ylabel('a6')

subplot(713)

plot(a5,'linewidth',2)

ylabel('a5')

subplot(714)

plot(a4,'linewidth',2)

ylabel('a4')

subplot(715)

plot(a3,'linewidth',2)

ylabel('a3')

subplot(716)

plot(a2,'linewidth',2)

ylabel('a2')

subplot(717)

plot(a1,'linewidth',2)

ylabel('a1')

xlabel('样本序号 N')%重构第1-7层细节系数

d7=wrcoef('d',c,l,'db5',7)

d6=wrcoef('d',c,l,'db5',6)

d5=wrcoef('d',c,l,'db5',5)

d4=wrcoef('d',c,l,'db5',4)

d3=wrcoef('d',c,l,'db5',3)

d2=wrcoef('d',c,l,'db5',2)

d1=wrcoef('d',c,l,'db5',1)

%显示细节系数

figure(4)

subplot(711)

plot(d7,'linewidth',2)

ylabel('d7')

subplot(712)

plot(d6,'linewidth',2)

ylabel('d6')

subplot(713)

plot(d5,'linewidth',2)

ylabel('d5')

subplot(714)

plot(d4,'linewidth',2)

ylabel('d4')

subplot(715)

plot(d3,'linewidth',2)

ylabel('d3')

subplot(716)

plot(d2,'linewidth',2)

ylabel('d2')

subplot(717)

plot(d1,'linewidth',2)

ylabel('d1')

xlabel('样本序号 N')

---