小波程序.txt

c++做的小波分析程序

用小波神经网络来对时间序列进行预测 

%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 variables 
clc; 
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 level 
      case 1, resolution = 1; 
      case 2, resolution = 2; 
      case 3, resolution = 3; 
      case 4, resolution = 4; 
end 

msg = ['Resolution level to be used is ', num2str(resolution)]; 
disp(msg); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% user desired amount of data to use 
data = 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 data 
      case 1, dataPoints = 48;        %1 day = 48 points 
      case 2, dataPoints = 96;        %2 days = 96 points 
      case 3, dataPoints = 144;       %3 days = 144 points 
      case 4, dataPoints = 192;       %4 days = 192 points 
      case 5, dataPoints = 240;       %5 days = 240 points 
      case 6, dataPoints = 288;       %6 days = 288 points 
      case 7, dataPoints = 336;       %1 weeks = 336 points 

end 

msg = ['No. of data points to be used is ', num2str(dataPoints)]; 
disp(msg); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%Menu for data set selection 
select = menu('Use QLD data of: ', 'Jan02',... 
                    'Feb02', 'Mar02 (Select this for testing)', 'Apr02', 'May02'); 
switch select 
      case 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 tDemandArray 
selectedDemandFile=[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 signal 
figure (1); 
subplot(resolution + 2, 1, 1); plot(tDemandArray(1: dataPoints)) 
legend('Demand original'); 
title('QLD Demand Data Signal'); 

%Plot the approximation demand data signal 
for i = 1 : resolution 
      subplot(resolution + 2, 1, i     + 1); plot(approx(1: dataPoints)) 
      legend('Demand Approximation'); 
end 

%After displaying approximation signal, display detail x 
for i = 1: resolution 
    if( i > 1 ) 
            detail(i, :) = dD(i-1, :)- dD(i, :); 
      else 
            detail(i, :) = tDemandArray' - dD(1, :); 
      end 

           if i == 1 
                 subplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints)) 
                 legendName = ['Demand Detail ', num2str(i)]; 
                  legend(legendName); 

            else 
                 subplot(resolution + 2, 1, resolution - i + 3); plot(detail(i, 1: dataPoints)) 
                 legendName = ['Demand Detail ', num2str(i)]; 
                 legend(legendName); 

            end 

    i = i + 1; 
end 

%Normalising approximation demand data 
maxDemand = max(approx'); %Find largest component 
normDemand = approx ./ maxDemand; %Right divison 

maxDemandDetail = [ ]; 
normDemandDetail = [, ]; 

detail = detail + 4000; 

for i = 1: resolution 
      maxDemandDetail(i) = max(detail(i, :)); 
      normDemandDetail(i, :) = detail(i, :) ./maxDemandDetail(i); 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%Reading the historical historical PRICE data into rrpArray 
selectedPriceFile = [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 signal 
figure (2); 
subplot(resolution + 3, 1, 1); plot(rrpArray(2: dataPoints)) 
legend('Price original'); 
title('Price Data Signal'); 

%Plot the approximation Price data signal 
for i = 1 : resolution 
      subplot(resolution + 3, 1, i + 1); plot(approxP(2: dataPoints)) 
      legend('Price Approximation'); 
end 

%After displaying approximation signal, display detail x 
for i = 1: resolution 
     if( i > 1 ) 
            detailP(i, :) = dP(i-1, :)- dP(i, :); 
      else 
            detailP(i, :) = rrpArray' - dP(1, :); 
      end 

if 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); 

else 
      subplot(resolution + 3, 1, resolution - i + 3); plot(detailP(i, 2: dataPoints)) 
      legendName = ['Price Detail ', num2str(i)]; 
      legend(legendName); 

end 
     i = i + 1; 
end 

%Normalising approximation Price data 
maxPrice = max(approxP'); %Find largest component 
normPrice = approxP ./ maxPrice; %Right divison 

maxPriceDetail = [ ]; 
normPriceDetail = [, ]; 

detailP = detailP + 40; 

for i = 1: resolution 
      maxPriceDetail(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 session 

while (1) 

      choice = menu('Please select one of the following options: ',... 
                          'Create new Neural Networks',... 
                          'Perform FORECASTING of load demand', 'QUIT session...'); 
      switch choice 
           case 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 NNs 
      if(scheme == '2') 
           nRetrain; 
      end 

      %If scheme is 'f', call <nForecast> to load the existing NN model 
      if(scheme == '3') 
           nForecast; 
      end 

      %If scheme is 'e', verifies and quit session if 'yes' is selected else continue 
      if(scheme == '4') 
           button = questdlg('Quit session?', 'Exit Dialog','Yes','No','No'); 
           switch button 
                 case 'Yes', disp(' '); 
                                   disp('Session has ended!!'); 
                                   disp(' '); 
                                   break; 
                 case 'No', quit cancel; 
           end 
      end 
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%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 2 
clc; 
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 signals 
for x = 1: numOfTimes 

      %Clearing variables 
      clear targetDemand; 
      clear inputs; 
      clear output; 
      clc; 

      if(x == 1) 
            neuralNetwork = ['Neural network settings for approximation level ', 
     num2str(resolution)]; 
      else 
            neuralNetwork = ['Neural network settings for detail level ', num2str(x - 1)]; 
      end 
      disp(neuralNetwork); 
      disp(' '); 

      %Set no. of input nodes and hidden neurons for the 
      %respective demand and price coefficient signal 
      numOfInputs = 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 examples 
     trainingLength = dataPoints; 

      %Set target outputs of the training examples 
     if(x == 1) 
           targetDemand = normDemand(targetStartAt: 1 + trainingLength); 
     else 
           targetDemand = 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 < trainingLength 
           if(x == 1) 
                inputs(1, y + 1) = normDemand(y + 1); 
                inputs(2, y + 1) = normPrice(y + 1); 

           else 
                inputs(1, y + 1) = normDemandDetail(x - 1, y + 1); 
                inputs(2, y + 1) = normPriceDetail(x - 1, y + 1); 

           end 

           output(y + 1, :) = targetDemand(y + 1); 

           y = y + 1; 
     end 

     inputs = (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 options 
      options = zeros(1, 18); 
      options(1) = 1; %Provides display of error values 
      options(14) = cycle * ni * np; 

     %Training the neural network 

      %netopt(net, options, x, t, alg); 
      nn = netopt(nn, options, inputs, output, 'scg'); 

      %Save the neural network 
      filename = ['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...'); 
      else 
            disp(' '); 
            disp('Model is now ready to forecast!!'); 
            disp(' '); 
            disp('Press ENTER key to continue...'); 
      end 

      pause; 
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 variables 
clc; 
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 NNs 
for x = 1 : numOfTimes 
      %Loading NN 
      filename = ['nn', num2str(x)]; 
      clear nn 
      load(filename); 
      clear in; 
      numOfInputs = nn.nin; 
      y = 0; 
      %Preparing details to forecast 
      while y < pointsAhead 
           if(x == 1) 
                 propData(1, y + 1) = normDemand(y + 1); 
                 propData(2, y + 1) = normPrice(y + 1); 

           else 
                 propData(1, y + 1) = normDemandDetail(x - 1, y + 1); 
                 propData(2, y + 1) = normPriceDetail(x - 1, y + 1); 

           end 

           y = y + 1; 
      end 

      if(x == 1) 
           forecast = mlpfwd(nn, propData') * maxDemand; 
      else 
            forecast = mlpfwd(nn, propData') * maxDemandDetail(x - 1) - 4000; 
      end 
      forecastResult(x, :) = forecast'; 
end 

%For debugging 
% forecastResult 
actualDemand = tDemandArray(2: 1 + pointsAhead); 
predicted = sum(forecastResult, 1)'; 

%Calculating Mean Absolute Error 
AbsError = abs(actualDemand - predicted(1: pointsAhead)) ./ actualDemand; 
msg = ['Mean Absolute Error = ', num2str(mean(AbsError(1: pointsAhead))), ' !!']; 
disp(' '); 
disp(msg); 

%Plot actual time series against predicted result 
figure(3) 
actualWithPredicted(:, 1) = actualDemand;