谐波治理检测和分析的Matlab程序

%---------------------说明--------------------------------------- %---------------------说明--------------------------------------- echo off; clear; %--------产生采样点:4周期每周期128点采样 nT = 10;nS = 128; %周期数和周期采样点数 NP = nS * nT; %总共的采样点数 t=[0.0 : 0.02 / nS : 0.02 * nT - 0.02 / nS]; %--------产生三相电流 i_ba = 100*sin(314*t); %A相基波,50Hz,380V,0相位 i_ha = 100*sin(314*3*t)/3+100*sin(314*5*t)/5+100*sin(314*7*t)/7+100*sin(314*9*t)/9+100*sin(314*11*t)/11; i_a = i_ba+i_ha; i_bb = 100*sin(314*t-2*pi/3); %A相基波,50Hz,380V,120相位 i_hb = 100*sin(314*3*t-2*pi/3)/3+100*sin(314*5*t-2*pi/3)/5+100*sin(314*7*t-2*pi/3)/7+100*sin(314*9*t-2*pi/3)/9+100*sin(314*11*t-2*pi/3)/11; i_b = i_bb+i_hb; i_bc = 100*sin(314*t+2*pi/3); %A相基波,50Hz,380V,-120相位 i_hc = 100*sin(314*3*t+2*pi/3)/3+100*sin(314*5*t+2*pi/3)/5+100*sin(314*7*t+2*pi/3)/7+100*sin(314*9*t+2*pi/3)/9+100*sin(314*11*t+2*pi/3)/11; i_c = i_bc+i_hc; %--------FFT变换 i_FFTa = fft (i_a, size (i_a, 2)); i_FFTb = fft (i_b, size (i_b, 2)); i_FFTc = fft (i_c, size (i_c, 2)); figure (1); subplot (3, 3, 1); plot (t, i_a, 'k'); xlabel ('时间t'); ylabel ('幅度'); title ('A相时域波形'); subplot (3, 3, 2); plot (t, i_b, 'k'); xlabel ('时间t'); ylabel ('幅度'); title ('B相时域波形'); subplot (3, 3, 3); plot (t, i_c, 'k'); xlabel ('时间t'); ylabel ('幅度'); title ('C相时域波形'); subplot (3, 3, 4); plot (t, abs (i_FFTa), 'k'); xlabel ('频率w'); ylabel ('幅度'); title ('A相频域波形'); subplot (3, 3, 5); plot (t, abs (i_FFTb), 'k'); xlabel ('频率w'); ylabel ('幅度'); title ('B相频域波形'); subplot (3, 3, 6); plot (t, abs (i_FFTc), 'k'); xlabel ('频率w'); ylabel ('幅度'); title ('C相频域波形'); subplot (3, 3, 7); plot (t, angle (i_FFTa), 'k'); xlabel ('频率w'); ylabel ('相度'); title ('A相相域波形'); subplot (3, 3, 8); plot (t, angle (i_FFTb), 'k'); xlabel ('频率w'); ylabel ('相度'); title ('B相相域波形'); subplot (3, 3, 9); plot (t, angle (i_FFTc), 'k'); xlabel ('频率w'); ylabel ('相度'); title ('C相相域波形'); %--------信号加入高斯白噪声 i_Ra = randn (1, size (i_a, 2)) * 20; %产生强度为10高斯白噪声 i_Rb = randn (1, size (i_b, 2)) * 20; i_Rc = randn (1, size (i_c, 2)) * 20; i_aR = i_a + i_Ra; %将噪声添加到信号中 i_bR = i_b + i_Rb; i_cR = i_c + i_Rc; i_FFTa = fft (i_aR, size (i_aR, 2)); %对加入噪声后的数据FFT变换 i_FFTb = fft (i_bR, size (i_bR, 2)); i_FFTc = fft (i_cR, size (i_cR, 2)); figure (2); subplot (3, 3, 1); plot (t, i_aR, 'k'); xlabel ('时间t'); ylabel ('幅度'); title ('A相加入噪声时域波形'); subplot (3, 3, 2); plot (t, i_bR, 'k'); xlabel ('时间t'); ylabel ('幅度'); title ('B相加入噪声时域波形'); subplot (3, 3, 3); plot (t, i_cR, 'k'); xlabel ('时间t'); ylabel ('幅度'); title ('C相加入噪声时域波形'); subplot (3, 3, 4); plot (t, abs (i_FFTa), 'k'); xlabel ('频率w'); ylabel ('幅度'); title ('A相频域波形'); subplot (3, 3, 5); plot (t, abs (i_FFTb), 'k'); xlabel ('频率w'); ylabel ('幅度'); title ('B相频域波形'); subplot (3, 3, 6); plot (t, abs (i_FFTc), 'k'); xlabel ('频率w'); ylabel ('幅度'); title ('C相频域波形'); subplot (3, 3, 7); plot (t, angle (i_FFTa), 'k'); xlabel ('频率w'); ylabel ('相度'); title ('A相相域波形'); subplot (3, 3, 8); plot (t, angle (i_FFTb), 'k'); xlabel ('频率w'); ylabel ('相度'); title ('B相相域波形'); subplot (3, 3, 9); plot (t, angle (i_FFTc), 'k'); xlabel ('频率w'); ylabel ('相度'); title ('C相相域波形'); %--------小波消噪 i_aS = WaveSoftDenoise8 (i_aR); %Db24小波8层到50Hz软消噪 i_aZ = WaveSoftDenoise3 (i_aS); %Db24小波3层到800Hz置零消噪 i_bS = WaveSoftDenoise8 (i_bR); %Db24小波8层到50Hz软消噪 i_bZ = WaveSoftDenoise3 (i_bS); %Db24小波3层到800Hz置零消噪 i_cS = WaveSoftDenoise8 (i_cR); %Db24小波8层到50Hz软消噪 i_cZ = WaveSoftDenoise3 (i_cS); %Db24小波3层到800Hz置零消噪 %-----------绘制去除噪声后的FFT i_FFTa = fft (i_aZ, size (i_aZ, 2)); %对消去噪声后的数据FFT变换 i_FFTb = fft (i_bZ, size (i_bZ, 2)); i_FFTc = fft (i_cZ, size (i_cZ, 2)); figure (3); subplot (3, 3, 1); plot (t, i_aZ, 'k'); xlabel ('时间t'); ylabel ('幅度'); title ('A相去除噪声时域波形'); subplot (3, 3, 2); plot (t, i_bZ, 'k'); xlabel ('时间t'); ylabel ('幅度'); title ('B相去除噪声时域波形'); subplot (3, 3, 3); plot (t, i_cZ, 'k'); xlabel ('时间t'); ylabel ('幅度'); title ('C相去除噪声时域波形'); subplot (3, 3, 4); plot (t, abs (i_FFTa), 'k'); xlabel ('频率w'); ylabel ('幅度'); title ('A相频域波形'); subplot (3, 3, 5); plot (t, abs (i_FFTb), 'k'); xlabel ('频率w'); ylabel ('幅度'); title ('B相频域波形'); subplot (3, 3, 6); plot (t, abs (i_FFTc), 'k'); xlabel ('频率w'); ylabel ('幅度'); title ('C相频域波形'); subplot (3, 3, 7); plot (t, angle (i_FFTa), 'k'); xlabel ('频率w'); ylabel ('相度'); title ('A相相域波形'); subplot (3, 3, 8); plot (t, angle (i_FFTb), 'k'); xlabel ('频率w'); ylabel ('相度'); title ('B相相域波形'); subplot (3, 3, 9); plot (t, angle (i_FFTc), 'k'); xlabel ('频率w'); ylabel ('相度'); title ('C相相域波形'); %--------小波在奇异点的诊断 i_aBn = randn (1, 5) * 20; %加入异常数据段 i_bBn = randn (1, 5) * 20; i_cBn = randn (1, 5) * 20; i_aB = i_a; i_bB = i_b; i_cB = i_c; for n = 1 : 5 %加入异常信号段导数据 i_aB (1, 128 * 2 + n)= i_a (1, 128 * 2 + n) + i_aBn (1, n); i_aB (1, 128 * 6 + n)= i_a (1, 128 * 6 + n) + i_aBn (1, n); i_bB (1, 128 * 2 + n)= i_b (1, 128 * 2 + n) + i_bBn (1, n); i_bB (1, 128 * 6 + n)= i_b (1, 128 * 6 + n) + i_bBn (1, n); i_cB (1, 128 * 2 + n)= i_c (1, 128 * 2 + n) + i_cBn (1, n); i_cB (1, 128 * 6 + n)= i_c (1, 128 * 6 + n) + i_cBn (1, n); end %--------小波分解奇异点判定 [Ca, La] = wavedec (i_aB, 3, 'db24'); [Cb, Lb] = wavedec (i_bB, 3, 'db24'); [Cc, Lc] = wavedec (i_cB, 3, 'db24'); cAa3 = appcoef (Ca, La, 'db24', 3); %取第4层低频系数 cDa3 = detcoef (Ca, La, 3); %取各层高频系数 cDa2 = detcoef (Ca, La, 2); cDa1 = detcoef (Ca, La, 1); cAb3 = appcoef (Cb, Lb, 'db24', 3); %取第4层低频系数 cDb3 = detcoef (Cb, Lb, 3); %取各层高频系数 cDb2 = detcoef (Cb, Lb, 2); cDb1 = detcoef (Cb, Lb, 1); cAc3 = appcoef (Cc, Lc, 'db24', 3); %取第4层低频系数 cDc3 = detcoef (Cc, Lc, 3); %取各层高频系数 cDc2 = detcoef (Cc, Lc, 2); cDc1 = detcoef (Cc, Lc, 1); i_aZ = WaveSoftDenoise3 (i_aB); %Db24小波3层到800Hz置零消噪 i_bZ = WaveSoftDenoise3 (i_bB); i_cZ = WaveSoftDenoise3 (i_cB); %--------绘制小波分解后的波形 figure (4); subplot (5, 3, 1); plot (i_aB, 'k'); title ('A相带奇异点波形'); axis([0, 1500, -200, 200]); subplot (5, 3, 2); plot (i_bB, 'k'); title ('B相带奇异点波形'); axis([0, 1500, -200, 200]); subplot (5, 3, 3); plot (i_cB, 'k'); title ('C相带奇异点波形'); axis([0, 1500, -200, 200]); subplot (5, 3, 4); plot (i_aZ, 'k'); title ('A相第3层低频信号'); subplot (5, 3, 5); plot (i_bZ, 'k'); title ('B相第3层低频信号'); subplot (5, 3, 6); plot (i_cZ, 'k'); title ('C相第3层低频信号'); subplot (5, 3, 7); plot (cDa3, 'k'); title ('A相第3层高频信号'); subplot (5, 3, 8); plot (cDb3, 'k'); title ('B相第3层高频信号'); subplot (5, 3, 9); plot (cDc3, 'k'); title ('C相第3层高频信号'); subplot (5, 3, 10); plot (cDa2, 'k'); title ('A相第2层高频信号'); subplot (5, 3, 11); plot (cDb2, 'k'); title ('B相第2层高频信号'); subplot (5, 3, 12); plot (cDb2, 'k'); title ('C相第2层高频信号'); subplot (5, 3, 13); plot (cDa1, 'k'); title ('A相第1层高频信号'); subplot (5, 3, 14); plot (cDb1, 'k'); title ('B相第1层高频信号'); subplot (5, 3, 15); plot (cDc1, 'k'); title ('C相第1层高频信号');