急求【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('时间')

用FFT就可以啦。加窗的效果作用是降低副瓣，同时意义上也是将两头的权值减小，以尽可能逼近真实情况。需要注意的是FFT之后的数据显示和你的采样率有关。

如果是实数采样，采样率为FS，则最大可检测频率为FS/2

如果是复数采样，采样率为FS，则最大可检测频率为FS。

检测的频率间隔为FS/N

我来回答吧

1.首先建立信号模型并采样

%谐波分析仿真数据生成 CreatData.m

clc

fs=3000f0=50

N=1024

n=1:N

t=(n-1)/fs

m=13

Am=[380.4 0 2.035 0 3.156 0 1.042 0 4.291 0 16.13 0 2.016 ]

PH=[0 0 60 0 135 0 157.5 0 0 0 60 0 0]

x=zeros(1,N)

for k=1 : m

x=x+Am(k)*cos(2*pi*f0*k*t+PH(k)/180*pi)

end

figure(1)

plot(x)

2.进行FFT谐波分析

%直接用FFT进行谐波计算 Harmonic_Analysis.m

N=length(x)

f=50

m=13

n=1:N

Fn=zeros(9,3)% 谐波频率

An=zeros(9,3)% 幅值

Pn=zeros(9,3)% 初相角

for i=1 : m

Fn(i,1) = iFn(i,2) = f * i

An(i,1) = iAn(i,2) = 0

Pn(i,1) = iPn(i,2) = 0

end

fsN=fs/N

f0=50

%1-FFT

X=fft(x)

X=X(1:N/2)

Xabs=abs(X)

for i= 1 : m

[Amax,index]=TriFind(Xabs,floor((i*f0-15)/fsN),ceil((i*f0+15)/fsN))

if(index==-1)

Fn(i,3) = 0

An(i,3) = 0

Pn(i,3) = 0

else

if(Xabs(index-1) >Xabs(index+1))

a1 = Xabs(index-1) / Xabs(index)

r1 = 1/(1+a1)

k01 = index -1

else

a1 = Xabs(index) / Xabs(index+1)

r1 = 1/(1+a1)

k01 = index

end

Fn(i,3) = (k01+r1-1)*fs/N

An(i,3) = 2*pi*r1*Xabs(k01)/(N*sin(r1*pi))

Pn(i,3) = phase(X(k01))-pi*r1

Pn(i,3) = Pn(i,3)*180/pi

Pn(i,3) = mod(Pn(i,3),180)

end

end

Fn

An

Pn

figure(2)

stem(Xabs(1:N/2)*2/N)

dlmwrite('An.txt',An,'delimiter', '\t','precision', '%.10f','newline', 'pc')

dlmwrite('Fn.txt',Fn,'delimiter', '\t','precision', '%.10f','newline', 'pc')

dlmwrite('Pn.txt',Pn,'delimiter', '\t','precision', '%.10f','newline', 'pc')

3.完全手打的，先运行CreatData.m，再运行Harmonic_Analysis.m，应该可行，有问题再补充

4.少了个函数文件

新建个.m,命名为 TriFind，将以下内容粘贴

function [ Amax,index ] = TriFind(A,first,last)

%TRIFIND find x[k] if x[k-1]<x[k] and x[k]>x[k+1]

% return x[k] and the k

Amax=0

index=-1

for k=first+1:last-1

if A(k-1)<A(k) &&A(k)>A(k+1)

Amax=A(k)

index=k

end

end

---