电力系统谐波分析代码

程序清单
（1）模拟谐波的信号的波形及频谱
fs=3000;%采样频率
N1=256;%采样点数
N2=1024;
n1=0:N1-1;
t=n1/fs;
x1=sin(2 * pi * t * 50)+0.3 * sin(2 * pi * t * 250)+0.2 * sin(2 * pi * t * 500);
%信号绘图部分
figure(1)
subplot(2,1,1);%绘制时域波形
plot(t,x1);
title(‘模拟谐波信号x1的波形’);
xlabel(‘时间t’);
w1=2000 * pi * (0:N2-1)/N2;
y1=x1 * exp(-1* i * t’ * w1);
subplot(2,1,2);%绘制频域域波形
plot(w1/(2 * pi),abs(y1));
title(‘模拟谐波信号x1的频谱’);
xlabel(‘频率f’);

（2）基波对应的FIR低通滤波器损耗函数曲线、基波信号波形及频谱
%50HZ的1次谐波信号
%利用FIR设计低通滤波器
wp=1/30 * 2 * pi;
ws=0.05 * 2 * pi;
Bt=ws-wp;%计算过渡带宽度
N0=ceil(6.6 * pi/Bt);
N=N0+mod(N0+1,2);%确保长度N是奇数
wc=(wp+ws)/2/pi;%关于pi归一化
hn=fir1(N,wc,‘low’,hamming(N+1));

%滤波器绘图部分
N2=1024;
hw=fft(hn,N2);
n=0:N;
k=1:N2/2;
w=2 * (0:N2/2-1)/N2;
figure(2)
subplot(2,1,1)
stem(n,hn);%绘制相应的滤波器波形
axis([0,200,-0.1,0.1]);
xlabel(‘n’);ylabel(‘h(n)’);
title(‘滤波器h(n)波形’)
grid on
subplot(2,1,2);%绘制相应滤波器的损耗函数曲线
plot(w,20 * log10(abs(hw(k))));
xlabel(‘w/π’);
ylabel(‘20lg|H_g(w)’);
title(‘损耗函数曲线’)
grid on

%经过低通滤波器后的信号
figure(3)
x2=conv(x1,hn);%时域卷积
fs=3000;%采样频率
Y1=length(x2);%采样后信号长度
n=0:Y1-1;
t=n/fs;
subplot(2,1,1);%绘制时域波形
plot(t,x2);
xlabel(‘时间t’);
title(‘基波信号波形’);
axis([0.03,0.10,-2,2]);
y2=x2 * exp(-1 * i * t’ * w1);
subplot(2,1,2);%绘制频域波形
plot(w1/(2 * pi),abs(y2));
xlabel(‘频率f’);
axis([0,600,0,150]);
title(‘基波信号频谱’);

（3）5次谐波对应的哈明窗损耗函数曲线、5次谐波信号波形及频谱
%50Hz的5次谐波信号
fs=3000;%采样频率
N1=256;%采样点数
N2=1024;
n1=0:N1-1;
t=n1/fs;
x1=sin(2 * pi * t * 50)+0.3 * sin(2 * pi * t * 250)+0.2 * sin(2 * pi * t * 500);

%利用FIR设计带通滤波器
wlp=1/15 * 2 * pi;wls=0.05 * 2 * pi;
wup=0.1 * 2 * pi ;wus=2/15 * 2 * pi;%设计指标参数赋值
Bt=wus-wup;
N0=ceil(6.6 * pi/Bt);
wp=[(wls+wlp)/2/pi,(wus+wup)/2/pi];%设置理想带通截止频率
N=N0+mod(N0+1,2); %确保长度N是奇数
hn=fir1(N,wp,‘bandpass’,hamming(N+1));

%滤波器绘图部分
figure(1)
hw=fft(hn,N2);
n=0:N;
k=1:N2/2;
w=2 * (0:N2/2-1)/N2;
subplot(2,1,1); %绘制哈明窗波形
stem(n,hn,’.’);
axis([0,100,-0.2,0.2]);
xlabel(‘n’);ylabel(‘h(n)’);
title(‘滤波器h(n)波形’)
grid on
subplot(2,1,2);%绘制相应滤波器h（n）的损耗函数曲线
plot(w,20 * log10(abs(hw(k))));
xlabel(‘w/π’);
ylabel(‘20lg|H_g(w)|’);
title(‘损耗函数曲线’)
axis([0,1,-150,50]);
grid on

%经过带通滤波器后的信号
x2=conv(x1,hn);%时域卷积
fs=3000;%采样频率
Y1=length(x2);%采样后信号长度
n=0:Y1-1;
t=n/fs;
figure(2)
subplot(2,1,1);%绘制时域波形
plot(t,x2);
xlabel(‘时间t’);
title(‘5次谐波信号波形’);
axis([0.02,0.1,-0.5,0.5]);
w=2000 * pi * (0:N2-1)/N2;
y2=x2 * exp(-1 * i * t’ * w);
subplot(2,1,2);%绘制频域波形
plot(w/(2 * pi),abs(y2));
xlabel(‘频率f’);
axis([0,600,0,50]);
title(‘5次谐波信号频谱’);

（4）10次谐波对应的哈明窗损耗函数曲线、10次谐波信号波形及频谱
%50Hz的10次谐波信号
fs=3000;%采样频率
N1=256;%采样点数
N2=1024;
n1=0:N1-1;
t=n1/fs;
x1=sin(2 * pi * t * 50)+0.3 * sin(2 * pi * t * 250)+0.2 * sin(2 * pi * t * 500);

%利用FIR设计高通滤波器
wp=0.15 * 2 * pi;
ws=2/15 * 2 * pi;
Bt=wp-ws;
N0=ceil(6.6 * pi/Bt);
N=N0+mod(N0+1,2);%确保h(n)长度N是奇数
wc=(wp+ws)/2/pi;%关于pi归一化
hn=fir1(N-1,wc,‘high’,hamming(N));

%滤波器绘图部分
figure(1)
hw=fft(hn,N2);
n=0:N-1;
k=1:N2/2;
w=2 * (0:N2/2-1)/N2;
subplot(2,1,1); %绘制哈明窗波形
stem(n,hn,’.’);
axis([0,200,-0.5,1]);
xlabel(‘n’);
ylabel(‘h(n)’);
title(‘滤波器h(n)波形’)
grid on
subplot(2,1,2);%绘制h（n）的损耗函数曲线
plot(w,20 * log10(abs(hw(k))));
axis([0,1,-150,5]);
xlabel(‘w/π’);
ylabel(‘20lg|H_g(w)|’);
title(‘损耗函数曲线’)
grid on

%经过高通滤波器后的信号
x2=conv(x1,hn);%时域卷积
fs=3000;%采样频率
Y1=length(x2);%采样后信号长度
n=0:Y1-1;
t=n/fs;
figure(2)%滤波后的信号绘图部分
subplot(2,1,1);%绘制时域波形
plot(t,x2);
xlabel(‘时间t’);
title(‘10次谐波信号波形’);
axis([0.04,0.1,-0.3,0.3]);
w=2000 * pi * (0:N2-1)/N2;
y2=x2 * exp(-1 * i * t’ * w);
subplot(2,1,2);%绘制频域波形
plot(w/(2 * pi),abs(y2));
xlabel(‘频率f’);
axis([0,1000,0,40]);
title(‘10次谐波信号频谱’);