www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/waveletpacketsdemo.m
%% Wavelet Packets: Decomposing the Details % This example shows how wavelet packets differ from the discrete wavelet % transform (DWT). The example shows how the wavelet packet transform % results in equal-width subband filtering of signals as opposed to the % coarser octave band filtering found in the DWT. % % This makes wavelet packets an attracive alternative to the DWT in a % number of applications. Two examples presented here are time-frequency % analysis and signal classification. % % If you use an orthogonal wavelet with the wavelet packet transform, you % additionally end up with a partioning of the signal energy among the % equal-width subbands. % % The example focuses on the 1-D case, but many of the concepts extend % directly to the wavelet packet transform of images. %% Discrete Wavelet and Discrete Wavelet Packet Transforms % The following figure shows a DWT tree for a 1-D input signal. % % <<../wv_tree.png>> % % *Figure 1:* DWT Tree down to level 3 for a 1-D input signal. % $\tilde{G}(f)$ is the scaling (lowpass) analysis filter and % $\tilde{H}(f)$ represents the wavelet (highpass) analysis filter. The % labels at the bottom show the partition of the frequency axis [0,1/2] % into subbands. % % The figure shows that subsequent levels of the DWT operate only on the % outputs of the lowpass (scaling) filter. At each level, the DWT divides % the signal into octave bands. In the critically-sampled DWT, the outputs % of the bandpass filters are downsampled by two at each level. In the % undecimated discrete wavelet transform, the outputs are not dowsampled. % % Compare the DWT tree with the full wavelet packet tree. % % <<../modwpt_nat_order.png>> % % *Figure 2:* Full wavelet packet tree down to level 3. % % In the wavelet packet transform, the filtering operations are also % applied to the wavelet, or detail, coefficients. The result is that % wavelet packets provide a subband filtering of the input signal into % progressively finer equal-width intervals. At each level, $j$, the % frequency axis [0,1/2] is divided into $2^j$ subbands. The subbands in % hertz at level j are approximately % % $[ \frac{n \mathrm{Fs}}{2^{j+1}}, \frac{(n+1) \mathrm{Fs}}{2^{j+1}}) \quad n=0,1,\ldots 2^j-1$ % where Fs is the sampling frequency. % % % How good this bandpass approximation is depends on how % frequency-localized the filters are. For Fejer-Korovkin filters like % 'fk18' (18 coefficients), the approximation is quite good. For a filter % like the Haar ('haar'), the approximation is not accurate. % % % In the critically-sampled wavelet packet transform the outputs of the % bandpass filters are downsampled by two. In the undecimated wavelet % packet transform, the outputs are not downsampled. %% Time-Frequency Analysis with Wavelet Packets % Because wavelet packets divide the frequency axis into finer intervals % than the DWT, wavelet packets are superior at time-frequency analysis. % As an example, consider two intermittent sine waves with frequencies of % 150 and 200 Hz in additive noise. The data are sampled at 1 kHz. To % prevent the loss of time resolution inherent in the critically-sampled % wavelet packet transform, use the undecimated transform. dt = 0.001; t = 0:dt:1-dt; x = ... cos(2*pi*150*t).*(t>=0.2 & t<0.4)+sin(2*pi*200*t).*(t>0.6 & t<0.9); y = x+0.05*randn(size(t)); [wpt,~,F] = modwpt(x,'TimeAlign',true); contour(t,F.*(1/dt),abs(wpt).^2); grid on; xlabel('Time (secs)'); ylabel('Hz'); title('Time-Frequency Analysis -- Undecimated Wavelet Packet Transform'); %% % Note that the wavelet packet transform is able to separate the 150 and % 200 Hz components. This is not true of the DWT, because 150 and 200 Hz % fall in the same octave band. The octave bands for a level-four DWT are % (in Hz) % % * [0,31.25) % * [31.25,62.5) % * [62.5,125) % * [125,250) % * [250,500) % % Demonstrate that these two components are time-localized by the DWT but % not separated in frequency. mra = modwtmra(modwt(y,'fk18',4),'fk18'); J = 5:-1:1; freq = 1/2*(1000./2.^J); contour(t,freq,flipud(abs(mra).^2)); grid on; ylim([0 500]) xlabel('Time (secs)'); ylabel('Hz'); title('Time-Frequency Analysis -- Undecimated Wavelet Transform'); %% % Of course, the continuous wavelet transform (CWT) provides a higher % resolution time-frequency analysis than the wavelet packet transform, but % wavelet packets have the added benefit of being an orthogonal transform % (when using an orthogonal wavelet). An orthogonal transform means that % the energy in the signal is preserved and partitioned among the % coefficients as the next section demonstrates. %% Energy Preservation in the Wavelet Packet Transform % In addition to filtering a signal into equal-width subbands at each % level, the wavelet packet transform partitions the signal's energy among % the subbands. This is true for both the decimated and undecimated wavelet % packet transforms. To demonstrate this, load a dataset containing a % seismic recording of an earthquake. This data is similar to the time % series used in the signal classification example below. load kobe; plot(kobe) grid on; xlabel('Seconds'); title('Kobe Earthquake Data'); %% % Obtain both the decimated and undecimated wavelet packet transforms of % the data down to level 3. To ensure consistent results for the decimated % wavelet packet transform, the example sets the |dwtmode| to periodization % and returns it to your original setting at the end of the example. st = dwtmode('status','nodisplay'); dwtmode('per','nodisplay'); wptreeDecimated = wpdec(kobe,3,'fk18'); wptUndecimated = modwpt(kobe,3); %% % First, extract all the wavelet packet coefficients at level three for the % decimated transform. wpcfs = zeros(8,381); for kk = 7:14 wpcfs(kk-6,:) = wpcoef(wptreeDecimated,kk); end %% % Now compute the total energy in both the decimated and undecimated % wavelet packet level-three coefficients and compare to the energy in the % original signal. decimatedEnergy = sum(sum(abs(wpcfs).^2,2)) undecimatedEnergy = sum(sum(abs(wptUndecimated).^2,2)) signalEnergy = norm(kobe,2)^2 %% % The DWT shares this important property with the wavelet packet transform. wt = modwt(kobe,'fk18',3); undecimatedWTEnergy = sum(sum(abs(wt).^2,2)) %% % Because the wavelet packet transform at each level divides the frequency % axis into equal-width intervals and partitions the signal energy among % those intervals, you can often construct useful feature vectors for % classification just by retaining the relative energy in each wavelet % packet. The next example demonstrates this. %% Wavelet Packet Classification -- Earthquake or Explosion? % Seismic recordings pick up activity from many sources. It is important to % be able to classify this activity based on its origin. For example, % earthquakes and explosions may present similar time-domain signatures, % but it is very important to differentiate between these two events. The % dataset |EarthQuakeExplosionData| contains 16 recordings with 2048 % samples. The first 8 recordings (columns) are seismic recordings of % earthquakes, the last 8 recordings (columns) are seismic recordings of % explosions. Load the data and plot an earthquake and explosion recording % for comparison. The data is taken from the R package astsa Stoffer % (2014), which supplements Shumway and Stoffer (2011). The author has % kindly allowed its use in this example. % % Plot a time series from each group for comparison. load EarthQuakeExplosionData; subplot(2,1,1) plot(EarthQuakeExplosionData(:,3)); xlabel('Time'); title('Earthquake Recording'); grid on; subplot(2,1,2); plot(EarthQuakeExplosionData(:,9)); xlabel('Time'); title('Explosion Recording'); grid on; %% % Form a wavelet packet feature vector by decomposing each time series down % to level three using the 'fk6' wavelet with an undecimated wavelet packet % transform. This results in 8 subbands with an approximate width of 1/16 % cycles/sample. Use the relative energy in each subband to create a % feature vector. As an example, obtain a wavelet packet relative energy % feature vector for the first earthquake recording. [wpt,~,F,E,RE] = modwpt(EarthQuakeExplosionData(:,1),3,'fk6'); %% % RE contains the relative energy in each of the 8 subbands. If you sum all % the elements in RE, it is equal to 1. Note that this is a significant % reduction in data. A time series of length 2048 is reduced to a vector % with 8 elements, each element representing the relative energy in the % wavelet packet nodes at level 3. The helper function % |helperEarthQuakeExplosionClassifier| obtains the relative energies for % each of the 16 recordings at level 3 using the 'fk6' wavelet. This % results in a 16-by-8 matrix and uses these feature vectors as inputs to a % kmeans classifier. The classifier uses the Gap statistic criterion to % determine the optimal number of clusters for the feature vectors between % 1 and 6 and classifies each vector. Additionally, the classifier performs % the exact same classification on the undecimated wavelet transform % coefficients at level 3 obtained with the 'fk6' wavelet and power spectra % for each of the time series. The undecimated wavelet transform results in % a 16-by-4 matrix (3 wavelet subbands and 1 scaling subband). The power % spectra result in a 16-by-1025 matrix. You must have the Statistics and % Machine Learning Toolbox(TM) and the Signal Processing Toolbox(TM) to run % the classifier. Level = 3; [WavPacketCluster,WtCluster,PspecCluster] = ... helperEarthQuakeExplosionClassifier(EarthQuakeExplosionData,Level) dwtmode(st,'nodisplay'); %% % |WavPacketCluster| is the clustering output for the undecimated wavelet % packet feature vectors. |WtCluster| is the clustering output using the % undecimated DWT feature vectors, and |PspecCluster| is the output for the % clustering based on the power spectra. The wavelet packet classification % has identified two clusters (|OptimalK: 2|) as the optimal number. % Examine the results of the wavelet packet clustering. WavPacketCluster.OptimalY %% % You see that only two errors have been made. Of the eight earthquake % recordings, seven are classified together (group 2) and one is mistakenly % classified as belonging to group 1. The same is true of the 8 explosion % recordings -- 7 are classified together and 1 is misclassified. % Classification based on the level-three undecimated DWT and power spectra % return one cluster as the optimal solution. % % If you repeat the above with the level equal to 4, the undecimated % wavelet transform and wavelet packet transform perform identically with % both identifying two clusters as optimal and missclassifying one time % series (the same two time series) in each group. %% Conclusions % In this example, you learned how the wavelet packet transform differs % from the discrete wavelet transform. Specifically, the wavelet packet % tree also filters the wavelet (detail) coefficients, while the wavelet % transform only iterates on the scaling (approximation) coefficients. % % You learned that the wavelet packet transform shares the important % energy-preserving property of the DWT while providing superior frequency % resolution. In some applications, this superior frequency resolution % makes the wavelet packet transform an attractive alternative to the DWT. %% References % Wickerhauser, M.V. "Adapted wavelet analysis from theory to software." % IEEE Press, 1994. % % Percival, D.B. and Walden, A.T. "Wavelet methods for time series % analysis". Cambridge University Press, 2000. % % Shumway, R.H. and Stoffer, D.H. "Time series analysis and its % applications with R examples." Springer, 2011. % % Stoffer, D.H. "astsa: Applied Statistical Time Series Analysis." R % package version 1.3. http://CRAN.R-project.org/package=astsa , 2014. %% Appendix % The following helper function is used in this example: % % * <matlab:edit('helperEarthQuakeExplosionClassifier') helperEarthQuakeExplosionClassifier>