用MATLAB编写的 小波程序

%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 %%%%%%%%%%%%%%%