小波变换程序

小波滤波器构造和消噪程序（2个）1．重构% mallet_wavelet.m% 此函数用于研究Mallet算法及滤波器设计% 此函数仅用于消噪a=pi/8; %角度赋初值b=pi/8;%低通重构FIR滤波器h0(n)冲激响应赋值h0=cos(a)*cos(b);h1=sin(a)*cos(b);h2=-sin(a)*sin(b);h3=cos(a)*sin(b);low_construct=[h0,h1,h2,h3];L_fre=4; %滤波器长度low_decompose=low_construct(end:-1:1); %确定h0(-n)，低通分解滤波器for i_high=1:L_fre; %确定h1(n)=(-1)^n,高通重建滤波器if(mod(i_high,2)==0);coefficient=-1;elsecoefficient=1;endhigh_construct(1,i_high)=low_decompose(1,i_high)*coefficient;endhigh_decompose=high_construct(end:-1:1); %高通分解滤波器h1(-n)L_signal=100; %信号长度n=1:L_signal; %信号赋值f=10;t=0.001;y=10*cos(2*pi*50*n*t).*exp(-20*n*t);figure(1);plot(y);title('原信号');check1=sum(high_decompose); %h0(n)性质校验check2=sum(low_decompose);check3=norm(high_decompose);check4=norm(low_decompose);l_fre=conv(y,low_decompose); %卷积l_fre_down=dyaddown(l_fre); %抽取,得低频细节h_fre=conv(y,high_decompose);h_fre_down=dyaddown(h_fre); %信号高频细节figure(2);subplot(2,1,1)plot(l_fre_down);title('小波分解的低频系数');subplot(2,1,2);plot(h_fre_down);title('小波分解的高频系数');l_fre_pull=dyadup(l_fre_down); %0差值h_fre_pull=dyadup(h_fre_down);l_fre_denoise=conv(low_construct,l_fre_pull);h_fre_denoise=conv(high_construct,h_fre_pull);l_fre_keep=wkeep(l_fre_denoise,L_signal); %取结果的中心部分,消除卷积影响h_fre_keep=wkeep(h_fre_denoise,L_signal);sig_denoise=l_fre_keep+h_fre_keep; %信号重构compare=sig_denoise-y; %与原信号比较figure(3);subplot(3,1,1)plot(y);ylabel('y'); %原信号subplot(3,1,2);plot(sig_denoise);ylabel('sig\_denoise'); %重构信号subplot(3,1,3);plot(compare);ylabel('compare'); %原信号与消噪后信号的比较2．消噪% mallet_wavelet.m% 此函数用于研究Mallet算法及滤波器设计% 此函数用于消噪处理%角度赋值%此处赋值使滤波器系数恰为db9%分解的高频系数采用db9较好,即它的消失矩较大%分解的有用信号小波高频系数基本趋于零%对于噪声信号高频分解系数很大，便于阈值消噪处理[l,h]=wfilters('db10','d');low_construct=l;L_fre=20; %滤波器长度low_decompose=low_construct(end:-1:1); %确定h0(-n)，低通分解滤波器for i_high=1:L_fre; %确定h1(n)=(-1)^n,高通重建滤波器if(mod(i_high,2)==0);coefficient=-1;elsecoefficient=1;endhigh_construct(1,i_high)=low_decompose(1,i_high)*coefficient;endhigh_decompose=high_construct(end:-1:1); %高通分解滤波器h1(-n) L_signal=100; %信号长度n=1:L_signal; %原始信号赋值f=10;t=0.001;y=10*cos(2*pi*50*n*t).*exp(-30*n*t);zero1=zeros(1,60); %信号加噪声信号产生zero2=zeros(1,30);noise=[zero1,3*(randn(1,10)-0.5),zero2];y_noise=y+noise;figure(1);subplot(2,1,1);plot(y);title('原信号');subplot(2,1,2);plot(y_noise);title('受噪声污染的信号');check1=sum(high_decompose); %h0(n),性质校验check2=sum(low_decompose);check3=norm(high_decompose);check4=norm(low_decompose);l_fre=conv(y_noise,low_decompose); %卷积l_fre_down=dyaddown(l_fre); %抽取,得低频细节h_fre=conv(y_noise,high_decompose);h_fre_down=dyaddown(h_fre); %信号高频细节figure(2);subplot(2,1,1)plot(l_fre_down);title('小波分解的低频系数');subplot(2,1,2);plot(h_fre_down);title('小波分解的高频系数');% 消噪处理for i_decrease=31:44;if abs(h_fre_down(1,i_decrease))>=0.000001h_fre_down(1,i_decrease)=(10^-7);endendl_fre_pull=dyadup(l_fre_down); %0差值h_fre_pull=dyadup(h_fre_down);l_fre_denoise=conv(low_construct,l_fre_pull);h_fre_denoise=conv(high_construct,h_fre_pull);l_fre_keep=wkeep(l_fre_denoise,L_signal); %取结果的中心部分,消除卷积影响h_fre_keep=wkeep(h_fre_denoise,L_signal);sig_denoise=l_fre_keep+h_fre_keep; %消噪后信号重构%平滑处理for j=1:2for i=60:70;sig_denoise(i)=sig_denoise(i-2)+sig_denoise(i+2)/2;end;end;compare=sig_denoise-y; %与原信号比较figure(3);subplot(3,1,1)plot(y);ylabel('y'); %原信号subplot(3,1,2);plot(sig_denoise);ylabel('sig\_denoise'); %消噪后信号subplot(3,1,3);plot(compare);ylabel('compare'); %原信号与消噪后信号的比较小波谱分析mallat算法经典程序clc;clear;%% 1.正弦波定义f1=50; % 频率1f2=100; % 频率2fs=2*(f1+f2); % 采样频率Ts=1/fs; % 采样间隔N=120; % 采样点数n=1:N;y=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts); % 正弦波混合figure(1)plot(y);title('两个正弦信号')figure(2)stem(abs(fft(y)));title('两信号频谱')%% 2.小波滤波器谱分析h=wfilters('db30','l'); % 低通g=wfilters('db30','h'); % 高通h=[h,zeros(1,N-length(h))]; % 补零（圆周卷积，且增大分辨率变于观察）g=[g,zeros(1,N-length(g))]; % 补零（圆周卷积，且增大分辨率变于观察）figure(3);stem(abs(fft(h)));title('低通滤波器图')figure(4);stem(abs(fft(g)));title('高通滤波器图')%% 3.MALLET分解算法(圆周卷积的快速傅里叶变换实现)sig1=ifft(fft(y).*fft(h)); % 低通(低频分量)sig2=ifft(fft(y).*fft(g)); % 高通(高频分量)figure(5); % 信号图subplot(2,1,1)plot(real(sig1));title('分解信号1')subplot(2,1,2)plot(real(sig2));title('分解信号2')figure(6); % 频谱图subplot(2,1,1)stem(abs(fft(sig1)));title('分解信号1频谱')subplot(2,1,2)stem(abs(fft(sig2)));title('分解信号2频谱')%% 4.MALLET重构算法sig1=dyaddown(sig1); % 2抽取sig2=dyaddown(sig2); % 2抽取sig1=dyadup(sig1); % 2插值sig2=dyadup(sig2); % 2插值sig1=sig1(1,[1:N]); % 去掉最后一个零sig2=sig2(1,[1:N]); % 去掉最后一个零hr=h(end:-1:1); % 重构低通gr=g(end:-1:1); % 重构高通hr=circshift(hr',1)'; % 位置调整圆周右移一位gr=circshift(gr',1)'; % 位置调整圆周右移一位sig1=ifft(fft(hr).*fft(sig1)); % 低频sig2=ifft(fft(gr).*fft(sig2)); % 高频sig=sig1+sig2; % 源信号%% 5.比较figure(7);subplot(2,1,1)plot(real(sig1));title('重构低频信号');subplot(2,1,2)plot(real(sig2));title('重构高频信号');figure(8);subplot(2,1,1)stem(abs(fft(sig1)));title('重构低频信号频谱');subplot(2,1,2)stem(abs(fft(sig2)));title('重构高频信号频谱');figure(9)plot(real(sig),'r','linewidth',2);hold on;plot(y);legend('重构信号','原始信号')title('重构信号与原始信号比较')使用小波包变换分析信号的MATLAB程序%t=0.001:0.001:1;t=1:1000;s1=sin(2*pi*50*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t));for t=1:500;s2(t)=sin(2*pi*50*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t));endfor t=501:1000;s2(t)=sin(2*pi*200*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t));endsubplot(9,2,1)plot(s1)title('原始信号')ylabel('S1')subplot(9,2,2)plot(s2)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]);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(s130);ylabel('S130');subplot(9,2,5);plot(s131);ylabel('S131');subplot(9,2,7);plot(s132);ylabel('S132');subplot(9,2,9);plot(s133);ylabel('S133');subplot(9,2,11);plot(s134);ylabel('S134');subplot(9,2,13);plot(s135);ylabel('S135');subplot(9,2,15);plot(s136);ylabel('S136');subplot(9,2,17);plot(s137);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]);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(s230);ylabel('S230');subplot(9,2,6);plot(s231);ylabel('S231');subplot(9,2,8);plot(s232);ylabel('S232');subplot(9,2,10);plot(s233);ylabel('S233');subplot(9,2,12);plot(s234);ylabel('S234');subplot(9,2,14);plot(s235);ylabel('S235');subplot(9,2,16);plot(s236);ylabel('S236');subplot(9,2,18);plot(s237);ylabel('S237');%fftfigurey1=fft(s1,1024);py1=y1.*conj(y1)/1024;y2=fft(s2,1024);py2=y2.*conj(y2)/1024;y130=fft(s130,1024);py130=y130.*conj(y130)/1024; y131=fft(s131,1024);py131=y131.*conj(y131)/1024; y132=fft(s132,1024);py132=y132.*conj(y132)/1024; y133=fft(s133,1024);py133=y133.*conj(y133)/1024; y134=fft(s134,1024);py134=y134.*conj(y134)/1024; y135=fft(s135,1024);py135=y135.*conj(y135)/1024; y136=fft(s136,1024);py136=y136.*conj(y136)/1024; y137=fft(s137,1024);py137=y137.*conj(y137)/1024; y230=fft(s230,1024);py230=y230.*conj(y230)/1024; y231=fft(s231,1024);py231=y231.*conj(y231)/1024; y232=fft(s232,1024);py232=y232.*conj(y232)/1024; y233=fft(s233,1024);py233=y233.*conj(y233)/1024; y234=fft(s234,1024);py234=y234.*conj(y234)/1024; y235=fft(s235,1024);py235=y235.*conj(y235)/1024; y236=fft(s236,1024);py236=y236.*conj(y236)/1024; y237=fft(s237,1024);py237=y237.*conj(y237)/1024; f=1000*(0:511)/1024; subplot(1,2,1);plot(f,py1(1:512));ylabel('P1');title('原始信号的功率谱') subplot(1,2,2);plot(f,py2(1:512));ylabel('P2');title('故障信号的功率谱') figuresubplot(4,2,1);plot(f,py130(1:512)); ylabel('P130');title('S130的功率谱') subplot(4,2,2);plot(f,py131(1:512)); ylabel('P131');title('S131的功率谱') subplot(4,2,3);plot(f,py132(1:512)); ylabel('P132'); subplot(4,2,4);plot(f,py133(1:512)); ylabel('P133'); subplot(4,2,5);plot(f,py134(1:512)); ylabel('P134'); subplot(4,2,6);plot(f,py135(1:512)); ylabel('P135'); subplot(4,2,7);plot(f,py136(1:512)); ylabel('P136'); subplot(4,2,8);plot(f,py137(1:512)); ylabel('P137'); figuresubplot(4,2,1);plot(f,py230(1:512)); ylabel('P230');title('S230的功率谱') subplot(4,2,2);plot(f,py231(1:512)); ylabel('P231');title('S231的功率谱') subplot(4,2,3);plot(f,py232(1:512)); ylabel('P232'); subplot(4,2,4);plot(f,py233(1:512)); ylabel('P233'); subplot(4,2,5);plot(f,py234(1:512)); ylabel('P234');subplot(4,2,6);plot(f,py235(1:512));ylabel('P235');subplot(4,2,7);plot(f,py236(1:512));ylabel('P236');subplot(4,2,8);plot(f,py237(1:512));ylabel('P237');figure%plottree(wpt)基于小波变换的图象去噪Normalshrink算法function [T_img,Sub_T]=threshold_2_N(img,levels)% reference :image denoising using wavelet thresholding[xx,yy]=size(img);HH=img((xx/2+1):xx,(yy/2+1):yy);delt_2=(std(HH(:)))^2;%(median(abs(HH(:)))/0.6745)^2;%T_img=img;for i=1:levelstemp_x=xx/2^i;temp_y=yy/2^i;% belt=1.0*(log(temp_x/(2*levels)))^0.5;belt=1.0*(log(temp_x/(2*levels)))^0.5; %2.5 0.8%HLHL=img(1:temp_x,(temp_y+1):2*temp_y);delt_y=std(HL(:));T_1=belt*delt_2/delt_y;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T_HL=sign(HL).*max(0,abs(HL)-T_1);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% T_img(1:temp_x,(temp_y+1):2*temp_y)=T_HL;Sub_T(3*(i-1)+1)=T_1;%LHLH=img((temp_x+1):2*temp_x,1:temp_y);delt_y=std(LH(:));T_2=belt*delt_2/delt_y;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%T_LH=sign(LH).*max(0,abs(LH)-T_2);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%T_img((temp_x+1):2*temp_x,1:temp_y)=T_LH;Sub_T(3*(i-1)+2)=T_2;%HHHH=img((temp_x+1):2*temp_x,(temp_y+1):2*temp_y);delt_y=std(HH(:));T_3=belt*delt_2/delt_y;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%T_HH=sign(HH).*max(0,abs(HH)-T_3);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %T_img((temp_x+1):2*temp_x,(temp_y+1):2*temp_y)=T_HH;Sub_T(3*(i-1)+3)=T_3;end用小波神经网络来对时间序列进行预测%File name : nprogram.m%Description : This file reads the data from its source into their respective matrices prior to% performing wavelet decomposition.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Clear command screen and variablesclc;clear;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % user desired resolution level (Tested: resolution = 2 is best)level = menu('Enter desired resolution level: ', '1',...'2 (Select this for testing)', '3', '4');switch levelcase 1, resolution = 1;case 2, resolution = 2;case 3, resolution = 3;case 4, resolution = 4;endmsg = ['Resolution level to be used is ', num2str(resolution)];disp(msg);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % user desired amount of data to usedata = menu('Choose amount of data to use: ', '1 day', '2 days', '3 days', '4 days',...'5 days', '6 days', '1 week (Select this for testing)');switch datacase 1, dataPoints = 48; %1 day = 48 pointscase 2, dataPoints = 96; %2 days = 96 pointscase 3, dataPoints = 144; %3 days = 144 pointscase 4, dataPoints = 192; %4 days = 192 pointscase 5, dataPoints = 240; %5 days = 240 pointscase 6, dataPoints = 288; %6 days = 288 pointscase 7, dataPoints = 336; %1 weeks = 336 pointsendmsg = ['No. of data points to be used is ', num2str(dataPoints)];disp(msg);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Menu for data set selectionselect = menu('Use QLD data of: ', 'Jan02',...'Feb02', 'Mar02 (Select this for testing)', 'Apr02', 'May02'); switch selectcase 1, demandFile = 'DATA200601_QLD1';case 2, demandFile = 'DATA200602_QLD1';case 3, demandFile = 'DATA200603_QLD1';case 4, demandFile = 'DATA200604_QLD1';case 5, demandFile = 'DATA200605_QLD1';end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Reading the historical DEMAND data into tDemandArrayselectedDemandFile=[demandFile,'.csv'];[regionArray, sDateArray, tDemandArray, rrpArray, pTypeArray] ...= textread(selectedDemandFile, '%s %q %f %f %s', 'headerlines', 1, 'delimiter', ',');%Display no. of points in the selected time series demand data [demandDataPoints, y] = size(tDemandArray);msg = ['The no. of points in the selected Demand data is ', num2str(demandDataPoints)]; disp(msg);%Decompose historical demand data signal[dD, l] = swtmat(tDemandArray, resolution, 'db2');approx = dD(resolution, :);%Plot the original demand data signalfigure (1);subplot(resolution + 2, 1, 1); plot(tDemandArray(1: dataPoints))legend('Demand Original');title('QLD Demand Data Signal');%Plot the approximation demand data signalfor i = 1 : resolutionsubplot(resolution + 2, 1, i + 1); plot(approx(1: dataPoints))legend('Demand Approximation');end%After displaying approximation signal, display detail xfor i = 1: resolutionif( i > 1 )detail(i, :) = dD(i-1, :)- dD(i, :);elsedetail(i, :) = tDemandArray' - dD(1, :);endif i == 1subplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints))legendName = ['Demand Detail ', num2str(i)];legend(legendName);elsesubplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints))legendName = ['Demand Detail ', num2str(i)];legend(legendName);endi = i + 1;end%Normalising approximation demand datamaxDemand = max(approx'); %Find largest componentnormDemand = approx ./ maxDemand; %Right divisonmaxDemandDetail = [ ];normDemandDetail = [, ];detail = detail + 4000;for i = 1: resolutionmaxDemandDetail(i) = max(detail(i, :));normDemandDetail(i, :) = detail(i, :) ./maxDemandDetail(i);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Reading the historical historical PRICE data into rrpArrayselectedPriceFile = [demandFile, '.csv'];[regionArray, sDateArray, tDemandArray, rrpArray, pTypeArray] ...= textread(selectedDemandFile, '%s %q %f %f %s', 'headerlines', 1, 'delimiter', ',');%Display no. of points in Price data[noOfDataPoints, y] = size(rrpArray);msg = ['The no. of points in Price data data is ', num2str(noOfDataPoints)];disp(msg);%Decompose historical Price data signal[dP, l] = swtmat(rrpArray, resolution, 'db2');approxP = dP(resolution, :);%Plot the original Price data signalfigure (2);subplot(resolution + 3, 1, 1); plot(rrpArray(2: dataPoints))legend('Price Original');title('Price Data Signal');%Plot the approximation Price data signalfor i = 1 : resolutionsubplot(resolution + 3, 1, i + 1); plot(approxP(2: dataPoints))legend('Price Approximation');end%After displaying approximation signal, display detail xfor i = 1: resolutionif( i > 1 )detailP(i, :) = dP(i-1, :)- dP(i, :);elsedetailP(i, :) = rrpArray' - dP(1, :);endif i == 1[B,A]=butter(1,0.65,'low');result =filter(B,A, detailP(i, 1: dataPoints));subplot(resolution + 3, 1, resolution - i + 4);plot(result(i, 2: dataPoints))legendName = ['low pass filter', num2str(i)];legend(legendName);subplot(resolution + 3, 1, resolution - i + 3); plot(detailP(i, 2: dataPoints)) legendName = ['Price Detail ', num2str(i)];legend(legendName);elsesubplot(resolution + 3, 1, resolution - i + 3); plot(detailP(i, 2: dataPoints)) legendName = ['Price Detail ', num2str(i)];legend(legendName);endi = i + 1;end%Normalising approximation Price datamaxPrice = max(approxP'); %Find largest componentnormPrice = approxP ./ maxPrice; %Right divisonmaxPriceDetail = [ ];normPriceDetail = [, ];detailP = detailP + 40;for i = 1: resolutionmaxPriceDetail(i) = max(detailP(i, :));normPriceDetail(i, :) = detailP(i, :) ./maxPriceDetail(i);end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Function here allows repetitive options to,% 1) Create new NNs, 2) Retrain the existing NNs,% 3) Perform load demand forecasting and 4) Quit sessionwhile (1)choice = menu('Please select one of the following options: ',...'CREATE new Neural Networks',...'Perform FORECASTING of load demand', 'QUIT session...');switch choicecase 1, scheme = '1';case 2, scheme = '2';case 3, scheme = '3';case 4, scheme = '4';end%If scheme is 'c', call <nCreate> to create new NNs, train them then perform forecast if(scheme == '1')nCreate;end%If scheme is 'r', call <nRetrain> to retrain the existing NNsif(scheme == '2')nRetrain;end%If scheme is 'f', call <nForecast> to load the existing NN modelif(scheme == '3')nForecast;end%If scheme is 'e', verifies and quit session if 'yes' is selected else continueif(scheme == '4')button = questdlg('Quit session?', 'Exit Dialog','Yes','No','No');switch buttoncase 'Yes', disp(' ');disp('Session has ended!!');disp(' ');break;case 'No', quit cancel;endendend%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%File name : ncreate.m%Description : This file prepares the input & output data for the NNs. It creates then trains the % NNs to mature them.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Clear command screen and set target demand ouput to start at point 2clc;targetStartAt = 2;disp('Program will now CREATE a Neural Network for training and forecasting...');disp(' ');disp('To capture the pattern of the signal, the model is ')disp('set to accept dataPoints x 2 sets of training examples; ');disp('1 set of demand + 1 sets of price. ');disp(' ');disp('The normalised demand data <point 2>, is to be taken as the ')disp('output value for the first iteration of training examples. ');disp(' ');disp('Press ENTER key to continue...');pause;%Clear command screen then prompt for no. of training cycles%For current program, 1 cycle = user set no. of iterations (ie: dataPoints)clc;cycle = input('Input the number of training cycles: ');numOfTimes = resolution + 1;%Creating and training the NNs for the respective%demand and price coefficient signalsfor x = 1: numOfTimes%Clearing variablesclear targetDemand;clear inputs;clear output;clc;if(x == 1)neuralNetwork = ['Neural network settings for approximation level ',num2str(resolution)];elseneuralNetwork = ['Neural network settings for detail level ', num2str(x - 1)];enddisp(neuralNetwork);disp(' ');%Set no. of input nodes and hidden neurons for the%respective demand and price coefficient signalnumOfInputs = 2;inputValue = ['Number of neural network INPUT units is set at ', num2str(numOfInputs)]; disp(inputValue);disp(' ');numOfOutput = 1;outValue = ['Output is set to ', num2str(numOfOutput)];disp(outValue);disp(' ');numOfHiddens = input('Enter the no. of HIDDEN units for the NN hidden : '); hiddenValue = ['Number of neural network HIDDEN units is set at ',num2str(numOfHiddens)];disp(hiddenValue);disp(' ');%Setting no. of training examplestrainingLength = dataPoints;%Set target outputs of the training examplesif(x == 1)targetDemand = normDemand(targetStartAt: 1 + trainingLength);elsetargetDemand = normDemandDetail(x - 1, targetStartAt: 1 + trainingLength);end%Preparing training examples%Setting training i/ps to be 2 with user defined no. of iterations (dataPoints)y = 0;while y < trainingLengthif(x == 1)inputs(1, y + 1) = normDemand(y + 1);inputs(2, y + 1) = normPrice(y + 1);elseinputs(1, y + 1) = normDemandDetail(x - 1, y + 1);inputs(2, y + 1) = normPriceDetail(x - 1, y + 1);endoutput(y + 1, :) = targetDemand(y + 1);y = y + 1;endinputs = (inputs');%Setting no. of training cycles[ni, np] = size(targetDemand); % <== [ni, np] tells the NN how long is 1 cycle; size(targetDemand)clear nn;%Create new neural network for respective coefficient signal%NET = MLP(NIN, NHIDDEN, NOUT, FUNC)nn = mlp(numOfInputs, numOfHiddens, numOfOutput, 'linear');%NN optionsoptions = zeros(1, 18);options(1) = 1; %Provides display of error valuesoptions(14) = cycle * ni * np;%Training the neural network%netopt(net, options, x, t, alg);nn = netopt(nn, options, inputs, output, 'scg');%Save the neural networkfilename = ['nn', num2str(x)];save(filename, 'nn');disp(' ');msg = ['Neural network successfully CREATED and saved as => ', filename]; disp(msg);if(x < 3)disp(' ');disp('Press ENTER key to continue training the next NN...');elsedisp(' ');disp('Model is now ready to forecast!!');disp(' ');disp('Press ENTER key to continue...');endpause;end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%File name : nforecast.m%Description : This file loads the trained NNs for load forcasting and %recombines the predicted % outputs from the NNs to form the final predicted load series.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Clear command screen and variablesclc;clear forecastResult;clear actualDemand;clear predicted;clear actualWithPredicted;disp('Neural networks loaded, performing electricity demand forecast...');disp(' ');pointsAhead = input('Enter no. of points to predict (should be < 336): ');numOfTimes = resolution + 1;%Predict coefficient signals using respective NNsfor x = 1 : numOfTimes%Loading NNfilename = ['nn', num2str(x)];clear nnload(filename);clear in;numOfInputs = nn.nin;y = 0;%Preparing details to forecastwhile y < pointsAheadif(x == 1)propData(1, y + 1) = normDemand(y + 1);propData(2, y + 1) = normPrice(y + 1);elsepropData(1, y + 1) = normDemandDetail(x - 1, y + 1);propData(2, y + 1) = normPriceDetail(x - 1, y + 1);endy = y + 1;endif(x == 1)forecast = mlpfwd(nn, propData') * maxDemand;elseforecast = mlpfwd(nn, propData') * maxDemandDetail(x - 1) - 4000;endforecastResult(x, :) = forecast';end%For debugging% forecastResultactualDemand = tDemandArray(2: 1 + pointsAhead);predicted = sum(forecastResult, 1)';%Calculating Mean Absolute ErrorAbsError = abs(actualDemand - predicted(1: pointsAhead)) ./ actualDemand;msg = ['Mean Absolute Error = ', num2str(mean(AbsError(1: pointsAhead))), ' !!']; disp(' ');disp(msg);%Plot actual time series against predicted resultfigure(3)actualWithPredicted(:, 1) = actualDemand;actualWithPredicted(:, 2) = predicted(1: pointsAhead);plot(actualWithPredicted);graph = ['Mean Absolute Error = ', num2str(mean(AbsError))];title(graph);legend('Actual', 'Forecasted'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %File name : nretrain.m%Description : This file loads the existing NNs and trains them again.%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clc;%Prompt for the starting point for trainingdisp('Program will now RETRAIN the Neural Networks ')disp('with the SAME intial data series again...');disp(' ');disp('To capture the pattern of the signal, the model is ')disp('set to accept dataPoints x 2 sets of training examples; ');disp('1 set of demand + 1 sets of price. ');disp(' ');disp('The normalised demand data <point 2>, is to be taken as the ')disp('output value for the first iteration of training examples. ');disp(' ');msg = ['Data points to be used for reTraining the NNs is from 1 to ',...num2str(dataPoints)];disp(msg);disp(' ');disp('Press ENTER key to continue...');pause;%Clear command screenclc;%Prompt for no. of training cycles%For current program, 1 cycle = user set no. of iterations (ie: dataPoints)cycle = input('Input number of cycles to retrain the NNs: ');numOfTimes = resolution + 1;%Loading existing NNs for trainingfor x = 1: numOfTimes%Re-initialising variablesclear targetDemand;clear inputs;clear output;clc;%Loading NN for the respective demand and temperature coefficient signals filename = ['nn', num2str(x)];clear nnload(filename);。