小波函数的选择如何影响cwt()的速度? [英] How does the choice of the wavelet function impact the speed of cwt()?

问题描述

在cwt()中，我可以指定要使用的小波函数.这对cwt()的速度有何影响?

In cwt() I can specify which wavelet function to use. How does that impact the speed of cwt()?

推荐答案

这里是一个基准，我使用 -singleCompThread 选项，当启动MATLAB强制它使用单个计算线程时. cwt()被传递了一个1,000,000个样本的信号，并被要求计算1到10的比例.我的CPU是i7-3610QM.

Here is a benchmark, which I run with the -singleCompThread option when starting MATLAB to force it to use a single computational thread. cwt() was passed a 1,000,000-sample signal and asked to compute scales 1 to 10. My CPU is an i7-3610QM.

使用的代码:

clear all

%% Benchmark parameters
results_file_name = 'results_scale1-10.csv';
number_of_random_runs = 10;
scales = 1:10;
number_of_random_samples = 1000000;

%% Construct a cell array containing all the wavelet names
wavelet_haar_names = {'haar'};
wavelet_db_names = {'db1'; 'db2'; 'db3'; 'db4'; 'db5'; 'db6'; 'db7'; 'db8'; 'db9'; 'db10'};
wavelet_sym_names = {'sym2'; 'sym3'; 'sym4'; 'sym5'; 'sym6'; 'sym7'; 'sym8'};
wavelet_coif_names = {'coif1'; 'coif2'; 'coif3'; 'coif4'; 'coif5'};
wavelet_bior_names = {'bior1.1'; 'bior1.3'; 'bior1.5'; 'bior2.2'; 'bior2.4'; 'bior2.6'; 'bior2.8'; 'bior3.1'; 'bior3.3'; 'bior3.5'; 'bior3.7'; 'bior3.9'; 'bior4.4'; 'bior5.5'; 'bior6.8'};
wavelet_rbior_names = {'rbio1.1'; 'rbio1.3'; 'rbio1.5'; 'rbio2.2'; 'rbio2.4'; 'rbio2.6'; 'rbio2.8'; 'rbio3.1'; 'rbio3.3'; 'rbio3.5'; 'rbio3.7'; 'rbio3.9'; 'rbio4.4'; 'rbio5.5'; 'rbio6.8'};
wavelet_meyer_names = {'meyr'};
wavelet_dmeyer_names = {'dmey'};
wavelet_gaus_names = {'gaus1'; 'gaus2'; 'gaus3'; 'gaus4'; 'gaus5'; 'gaus6'; 'gaus7'; 'gaus8'};
wavelet_mexh_names = {'mexh'};
wavelet_morl_names = {'morl'};
wavelet_cgau_names = {'cgau1'; 'cgau2'; 'cgau3'; 'cgau4'; 'cgau5'};
wavelet_shan_names = {'shan1-1.5'; 'shan1-1'; 'shan1-0.5'; 'shan1-0.1'; 'shan2-3'};
wavelet_fbsp_names = {'fbsp1-1-1.5'; 'fbsp1-1-1'; 'fbsp1-1-0.5'; 'fbsp2-1-1'; 'fbsp2-1-0.5'; 'fbsp2-1-0.1'};
wavelet_cmor_names = {'cmor1-1.5'; 'cmor1-1'; 'cmor1-0.5'; 'cmor1-1'; 'cmor1-0.5'; 'cmor1-0.1'};

% Concatenate all wavelet names into a single cell array 
wavelet_categories_names = who('wavelet*names');
wavelet_names = {};
for wavelet_categories_number=1:size(wavelet_categories_names,1)
    temp = wavelet_categories_names(wavelet_categories_number);
    temp = eval(temp{1});
    wavelet_names = vertcat(wavelet_names, temp);
end

%% Prepare data
random_signal = rand(number_of_random_runs,number_of_random_samples); 


%% Run benchmarks
result_file_ID = fopen(results_file_name, 'w');
for wavelet_number = 1:size(wavelet_names,1)
   wavelet_name = wavelet_names(wavelet_number,:)

   % Compute wavelet on a random signal
   tic
   for run = 1:number_of_random_runs             
       cwt(random_signal(run, :),scales,wavelet_name{1});       
   end
   run_time_random_test = toc

   fprintf(result_file_ID, '%s,', wavelet_name{1})
   fprintf(result_file_ID, '%d\n', run_time_random_test)
end

size(wavelet_names,1)
fclose(result_file_ID);

如果要查看比例尺的影响:

If you want to see the impact of the choice of the scale:

使用的代码:

clear all

%% Benchmark parameters
results_file_name = 'results_sym2_change_scale.csv';
number_of_random_runs = 10;
scales = 1:10;
number_of_random_samples = 10000000;
% wavelet_names = {'sym2', 'sym3'}%, 'sym4'};
output_directory = 'output';
wavelet_names = get_all_wavelet_names();

%% Prepare data
random_signal = rand(number_of_random_runs,number_of_random_samples);

%% Prepare result folder
if ~exist(output_directory, 'dir')
  mkdir(output_directory);
end

%% Run benchmarks
result_file_ID = fopen(results_file_name, 'w');
for wavelet_number = 1:size(wavelet_names,1)
    wavelet_name = wavelet_names{wavelet_number}
    if wavelet_number > 1
        fprintf(result_file_ID, '%s\n', '');
    end
    fprintf(result_file_ID, '%s', wavelet_name)
    run_time_random_test_scales = zeros(size(scales,2),1);
    for scale_number = 1:size(scales,2)
       scale = scales(scale_number);
       % Compute wavelet on a random signal
       tic
       for run = 1:number_of_random_runs             
           cwt(random_signal(run, :),scale,wavelet_name);       
       end
       run_time_random_test = toc
       fprintf(result_file_ID, ',%d', run_time_random_test)
       run_time_random_test_scales(scale_number) = run_time_random_test;
    end
    figure
    bar(run_time_random_test_scales)
    title(['Run time on random signal for ' wavelet_name])
    xlabel('Scale')
    ylabel('Run time (seconds)')
    save_figure( fullfile(output_directory, ['run_time_random_test_' wavelet_name]) )
    close all 
end

size(wavelet_names,1)
fclose(result_file_ID);

具有3个功能:

get_all_wavelet_names.m:

function [ wavelet_names ] = get_all_wavelet_names(  )
%GET_ALL_WAVELET_NAMES Get a list of available wavelet functions

%% Construct a cell array containing all the wavelet names
wavelet_haar_names = {'haar'};
wavelet_db_names = {'db1'; 'db2'; 'db3'; 'db4'; 'db5'; 'db6'; 'db7'; 'db8'; 'db9'; 'db10'};
wavelet_sym_names = {'sym2'; 'sym3'; 'sym4'; 'sym5'; 'sym6'; 'sym7'; 'sym8'};
wavelet_coif_names = {'coif1'; 'coif2'; 'coif3'; 'coif4'; 'coif5'};
wavelet_bior_names = {'bior1.1'; 'bior1.3'; 'bior1.5'; 'bior2.2'; 'bior2.4'; 'bior2.6'; 'bior2.8'; 'bior3.1'; 'bior3.3'; 'bior3.5'; 'bior3.7'; 'bior3.9'; 'bior4.4'; 'bior5.5'; 'bior6.8'};
wavelet_rbior_names = {'rbio1.1'; 'rbio1.3'; 'rbio1.5'; 'rbio2.2'; 'rbio2.4'; 'rbio2.6'; 'rbio2.8'; 'rbio3.1'; 'rbio3.3'; 'rbio3.5'; 'rbio3.7'; 'rbio3.9'; 'rbio4.4'; 'rbio5.5'; 'rbio6.8'};
wavelet_meyer_names = {'meyr'};
wavelet_dmeyer_names = {'dmey'};
wavelet_gaus_names = {'gaus1'; 'gaus2'; 'gaus3'; 'gaus4'; 'gaus5'; 'gaus6'; 'gaus7'; 'gaus8'};
wavelet_mexh_names = {'mexh'};
wavelet_morl_names = {'morl'};
wavelet_cgau_names = {'cgau1'; 'cgau2'; 'cgau3'; 'cgau4'; 'cgau5'};
wavelet_shan_names = {'shan1-1.5'; 'shan1-1'; 'shan1-0.5'; 'shan1-0.1'; 'shan2-3'};
wavelet_fbsp_names = {'fbsp1-1-1.5'; 'fbsp1-1-1'; 'fbsp1-1-0.5'; 'fbsp2-1-1'; 'fbsp2-1-0.5'; 'fbsp2-1-0.1'};
wavelet_cmor_names = {'cmor1-1.5'; 'cmor1-1'; 'cmor1-0.5'; 'cmor1-1'; 'cmor1-0.5'; 'cmor1-0.1'};

% Concatenate all wavelet names into a single cell array 
wavelet_categories_names = who('wavelet*names');
wavelet_names = {};
for wavelet_categories_number=1:size(wavelet_categories_names,1)
    temp = wavelet_categories_names(wavelet_categories_number);
    temp = eval(temp{1});
    wavelet_names = vertcat(wavelet_names, temp);
end


end

save_figure.m:

function [  ] = save_figure( output_graph_filename )
% Record aa figure as PNG and fig files

% Create the folder if it doesn't exist already.
[pathstr, name, ext] = fileparts(output_graph_filename);
if ~exist(pathstr, 'dir')
    mkdir(pathstr);
end

h = gcf;
set(0,'defaultAxesFontSize',18) % http://www.mathworks.com/support/solutions/en/data/1-8XOW94/index.html?solution=1-8XOW94
boldify(h);
print('-dpng','-r600', [output_graph_filename '.png']);
print(h,[output_graph_filename '.pdf'],'-dpdf','-r600')
saveas(gcf,[output_graph_filename '.fig'], 'fig')

end

和boldify.m:

function boldify(h,g)
%BOLDIFY        Make lines and text bold for standard viewgraph style.
%               BOLDIFY boldifies the lines and text of the current figure.
%               BOLDIFY(H) applies to the graphics handle H.
%
%               BOLDIFY(X,Y) specifies an X by Y inch graph of the current
%               figure.  If text labels have their 'UserData' data property
%               set to 'slope = ...', then the 'Rotation' property is set to
%               account for changes in the graph's aspect ratio.  The
%               default is MATLAB's default.
%               S. T. Smith
% The name of this function does not represent an endorsement by the author
% of the egregious grammatical trend of verbing nouns.

if nargin < 1, h = gcf;, end

% Set (and get) the default MATLAB paper size and position
set(gcf,'PaperPosition','default');
units = get(gcf,'PaperUnits');
set(gcf,'PaperUnits','inches');
fsize = get(gcf,'PaperPosition');
fsize = fsize(3:4);                     % Figure size (X" x Y") on paper.
psize = get(gcf,'PaperSize');

if nargin == 2                          % User specified graph size
 fsize = [h,g];
 h = gcf;
end

% Set the paper position of the current figure
set(gcf,'PaperPosition', ...
 [(psize(1)-fsize(1))/2 (psize(2)-fsize(2))/2 fsize(1) fsize(2)]);
fsize = get(gcf,'PaperPosition');
fsize = fsize(3:4);                     % Graph size (X" x Y") on paper.

set(gcf,'PaperUnits',units);            % Back to original

% Get the normalized axis position of the current axes
units = get(gca,'Units');
set(gca,'Units','normalized');
asize = get(gca,'Position');
asize = asize(3:4);
set(gca,'Units',units);
ha = get(h,'Children');

for i=1:length(ha)
%   if get(ha(i),'Type') == 'axes'
   % changed by B. A. Miller
 if strcmp(get(ha(i), 'Type'), 'axes') == 1
   units = get(ha(i),'Units');
   set(ha(i),'Units','normalized');
   asize = get(ha(i),'Position'); % Axes Position (normalized)
   asize = asize(3:4);
   set(ha(i),'Units',units);

   [m,j] = max(asize); j = j(1);
   scale = 1/(asize(j)*fsize(j));      % scale*inches -normalized units
   set(ha(i),'FontWeight','Bold');
   set(ha(i),'LineWidth',2);

   [m,k] = min(asize); k = k(1);
   if asize(k)*fsize(k) > 1/2
     set(ha(i),'TickLength',[1/8 1.5*1/8]*scale); % Gives 1/8" ticks
   else
     set(ha(i),'TickLength',[3/32 1.5*3/32]*scale); % Gives 3/32" ticks
   end

   set(get(ha(i),'XLabel'),'FontSize',18); % 14-pt labels
   set(get(ha(i),'XLabel'),'FontWeight','Bold');
   set(get(ha(i),'XLabel'),'VerticalAlignment','top');

   set(get(ha(i),'YLabel'),'FontSize',18); % 14-pt labels
   set(get(ha(i),'YLabel'),'FontWeight','Bold');
   %set(get(ha(i),'YLabel'),'VerticalAlignment','baseline');

   set(get(ha(i),'Title'),'FontSize',18); % 16-pt titles
   set(get(ha(i),'Title'),'FontWeight','Bold');
%    set(get(ha(i), 'FontSize',20, 'XTick',[]));
 end
 hc = get(ha(i),'Children');
 for j=1:length(hc)
   chtype = get(hc(j),'Type');
   if chtype(1:4) == 'text'
     set(hc(j),'FontSize',17);         % 12 pt descriptive labels
     set(hc(j),'FontWeight','Bold');
     ud = get(hc(j),'UserData');       % User data
     if length(ud) 8
       if ud(1:8) == 'slope = '        % Account for change in actual slope
         slope = sscanf(ud,'slope = %g');
         slope = slope*(fsize(2)/fsize(1))/(asize(2)/asize(1));
         set(hc(j),'Rotation',atan(slope)/pi*180);
       end
     end
   elseif chtype(1:4) == 'line'
     set(hc(j),'LineWidth',2);
   end
 end
end

---

奖金:前100个样本中有1000000个样本的随机信号上所有小波之间的相关性:

---

Bonus: correlation between all wavelets on a random signal with 1000000 samples with the first 10 scales:

使用的代码:

%% PRE-REQUISITE: You need to download http://www.mathworks.com/matlabcentral/fileexchange/24253-customizable-heat-maps , which gives the function heatmap()

%% Benchmark parameters
scales = 1:10;
number_of_random_samples = 1000000;
% wavelet_names = {'sym2'; 'sym3'; 'sym4'; 'sym5'; 'sym6'; 'sym7'; 'sym8'};
% wavelet_names = {'cgau1'; 'cgau2'; 'cgau3'; 'cgau4'; 'cgau5'};
wavelet_names = {'db2'; 'sym2'};

OUTPUT_FOLDER = 'output_corr';
% wavelet_names = get_all_wavelet_names(); % WARNING: you need to remove all complex wavelets, viz. cgau1, shan, fbsp and cmor, and the heatmap will be pissed to see complex values coming to her.

%% Prepare data
random_signal = rand(1,number_of_random_samples);
results = zeros(size(wavelet_names,1), number_of_random_samples);

%% Prepare result folder
if ~exist(OUTPUT_FOLDER, 'dir')
    mkdir(OUTPUT_FOLDER);
end

%% Run benchmarks
for scale_number = 1:size(scales,2)
    scale = scales(scale_number);
    for wavelet_number = 1:size(wavelet_names,1)
        wavelet_name = wavelet_names{wavelet_number}

        % Compute wavelet on a random signal
        run = 1;
        results(wavelet_number, :) = cwt(random_signal(run, :),scale,wavelet_name);

        if wavelet_number == 999
            break
        end
    end

    correlation_results = corrcoef(results')
    heatmap(correlation_results, [], [], '%0.2f', 'MinColorValue', -1.0, 'MaxColorValue', 1.0, 'Colormap', 'jet',...
        'Colorbar', true, 'ColorLevels', 64, 'UseFigureColormap', false);

    title(['Correlation matrix for scale ' num2str(scale)]);
    xlabel(['Wavelet 1 to ' num2str(size(wavelet_names,1)) ' for scale ' num2str(scale)]);
    ylabel(['Wavelet 1 to ' num2str(size(wavelet_names,1)) ' for scale ' num2str(scale)]);
    snapnow
    print('-dpng','-r600',fullfile(OUTPUT_FOLDER, ['scalecorr' num2str(scale) '.png']))

end

每个小波在不同尺度(1到100)之间的相关性:

Correlation for each wavelet between different scales (1 to 100):

使用的代码:

%% PRE-REQUISITE: You need to download http://www.mathworks.com/matlabcentral/fileexchange/24253-customizable-heat-maps , which gives the function heatmap()

%% Benchmark parameters
scales = 1:100;
number_of_random_samples = 1000000;
% wavelet_name = 'gaus2';
% wavelet_names = {'sym2', 'sym3'}%, 'sym4'};
OUTPUT_FOLDER = 'output_corr';
wavelet_names = get_all_wavelet_names();  % WARNING: you need to remove all complex wavelets, viz. cgau1, shan, fbsp and cmor, and the heatmap will be pissed to see complex values coming to her.

%% Prepare data
random_signal = rand(1,number_of_random_samples);
results = zeros(size(scales,2), number_of_random_samples); 

%% Prepare result folder
if ~exist(OUTPUT_FOLDER, 'dir')
  mkdir(OUTPUT_FOLDER);
end

%% Run benchmarks
for wavelet_number = 1:size(wavelet_names,1)
    wavelet_name = wavelet_names{wavelet_number}
    run_time_random_test_scales = zeros(size(scales,2),1);
    for scale_number = 1:size(scales,2)
       scale = scales(scale_number);
       run = 1;
       % Compute wavelet on a random signal
       results(scale_number, :) = cwt(random_signal(run, :),scale,wavelet_name);       
    end

   correlation_results = corrcoef(results')    
    heatmap(correlation_results, [], [], '%0.2f', 'MinColorValue', -1.0, 'MaxColorValue', 1.0, 'Colormap', 'jet',...
            'Colorbar', true, 'ColorLevels', 64, 'UseFigureColormap', false);

    title(['Correlation matrix for wavelet ' wavelet_name]);
    xlabel(['Scales 1 to ' num2str(max(scales)) ' for wavelet ' wavelet_name]);
    ylabel(['Scales 1 to ' num2str(max(scales)) ' for wavelet ' wavelet_name]);
    snapnow
    print('-dpng','-r600',fullfile(OUTPUT_FOLDER, [wavelet_name '_scalecorr_scale1to' num2str(max(scales)) '.png']))
end

这篇关于小波函数的选择如何影响cwt()的速度?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！