www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/HaarTransformsForTimeSeriesDataAndImagesExample.m
%% Haar Transforms for Time Series Data and Images % This example shows how to use Haar transforms to analyze time series % data and images. To run all the code in this example, you must have % the Signal Processing Toolbox(TM) and Image Processing Toolbox(TM). %% % First, visualize the Haar wavelet. [~,psi,x] = wavefun('haar',10); x = x(2:end-1); psi = psi(2:end-1); hl = plot(x(1:512),psi(1:512)); grid on; hold on; line(x(513:end),psi(513:end)); xlabel('t'); ylabel('\psi(t)','fontsize',14); title('Haar Wavelet'); %% % The Haar wavelet is discontinuous. As a result, it is typically not used % in denoising or compression applications where smoothness of the % reconstruction wavelet is an important consideration. However, Haar % transforms are useful in a number of applications due to their superior % time (spatial) localization and computational efficiency. The % Wavelet Toolbox(TM) supports Haar analysis in most of the discrete % wavelet analysis tools. This example features Haar lifting % implementations which support integer-to-integer wavelet transforms for % both 1-D and 2-D data and multichannel (multivariate) 1-D data. %% Analyze Signal Variability By Scale % Load and plot the |clock_571| dataset. This example is essentially a % recreation of the analysis presented in Percival & Walden (2000), pp % 13-16. load clock_571; figure; plot(clock_571) xlabel('Days') grid on; title('Daily Average Fractional Frequency Deviates -- Cesium Clock'); %% % The data are daily average fractional frequency deviates for a particular % cesium beam atomic clock with respect to the U.S. Naval Observatory % master clock. If the value of the time series is 0, that means the cesium % clock has neither lost nor gained time with respect to the master clock % over the duration of the day. If the value is negative, the clock has % lost time that day, a postive value means that the clock has gained time. % For this data, the values are all negative. For certain applications, % like geodesy, it is important to know whether there are certain time % scales where the clock's deviation from the master clock is at its lowest % value. In other words, are there certain scales where the clock agrees % most closely with the master clock? The Haar transform is useful here % because it possesses two important properties: It decorrelates data by % scale and it partitions signal energy among scale. %% % To illustrate the decorrelating property, obtain the Haar transform down % to level 6. Plot the autocorrelation sequence of the original data along % with the autocorrelation of the wavelet coefficients by scale for scales % of 2,4,8,16, and 32 days. Dashed lines on the plots delineate 95% % confidence intervals for white noise inputs. Values exceeding those lines % are indicative of significant autocorrelation in the data. [s,w] = haart(clock_571,6); helperAutoCorr(clock_571,w); %% % The top plot shows the autocorrelation sequence for the original data. % Subsequent plots show the autocorrelation sequences for wavelet % coefficients at increasingly coarser scales. It is clear that the % autocorrelation sequence of the original data exhibits correlation at all % lags while the Haar transform coefficients are decorrelated. Next, % demonstrate energy conservation. sigenergy = norm(clock_571,2)^2 energyByScale = cellfun(@(x)norm(x,2)^2,w); haarenergy = norm(s,2)^2+sum(energyByScale) %% % The total signal energy is preserved by the Haar transform. Because of % these properties, you can make meaningful inferences based on the % proportion of signal energy captured by the wavelet coefficients at each % scale. scales = 2.^(1:6); figure; plot(scales,energyByScale,'-o') xlabel('Scale (days)'); set(gca,'xscale','log') set(gca,'xtick',2.^(1:6)); ylabel('Proportion of Signal Energy'); grid on; %% % You see that the energy is at a minimum for scales of 16 and 32 days. For % the Haar wavelet (and all Daubechies wavelets), the wavelet coefficients % at a given scale represent differences between weighted averages of the % data over a duration 1/2 the length of the scale. This plot indicates % scales over which the cesium clock is in best agreement with the master % clock. This means that that considering data over approximately two week % or even one month periods is more accurate than data on smaller or longer % scales. As previously mentioned, this has important implications for % geodesy where extremely precise time measurements are critical. %% % Although the Haar wavelet is discontinuous, it is still effective at % representing various kinds of time series. Examples include count data % and data where values of a time series are tied to some specific state, % which affects the level of the time series. As an example, consider the % relationship between heart rate and sleep state. %% Create Signal Approximations % The data consist of two time series. One time series is the heart rate of % a 66-day old infant sampled every 16 seconds for just over 9 hours. The % heart rate time series is integer-valued. The other time series is the % expertly scored sleep state of the same infant over the same period with % the same sampling rate. % The sleep state data was scored based on the infant's EEG and % EOG (eye movement) data, not based on the heart rate. % The sleep state codes are 1=quiet sleep, 2=between quiet and active % sleep, 3=active sleep, and 4=awake. % Both time series were recorded by Prof. Peter Fleming, Dr Andrew % Sawczenko, and Jeanine Young of the Institute of Child Health, Royal % Hospital for Sick Children, Bristol and kindly provided for use in this % example. Plot the heart rate data along with the sleep states. load BabyECGData; figure; yyaxis left; plot(times,HR); ylabel('HR'); xlabel('Hrs'); YLim = [min(HR)-1 max(HR)+1]; yyaxis right; plot(times,SS); ylabel('Sleep State'); YLim = [0.5 4.5]; title('Baby ECG and Sleep State'); %% % An inspection of the data reveals an apparent correlation between sleep % state and heart rate, but the data is extremely noisy. Because the Haar % transform provides a staircase approximation to a signal, it is often % useful in situations where a response is dependent on a predictor % variable with a small number of discrete states. Here the discrete states % are the four sleep stages. Obtain the Haar approximation of the heart % rate data using a level 5 approximation. Because the heart rate data is % integer-valued, use the |'integer'| flag to ensure that integer-valued % data is returned. Plot the result. [S,W] = haart(HR,'integer'); HaarHR = ihaart(S,W,5,'integer'); figure; hL = plotyy(times,HaarHR,times,SS); Ax1 = hL(1); Ax2 = hL(2); Ax1.YLim = [min(HaarHR)-1 max(HaarHR)+1]; Ax1.YLabel.String = 'HR'; Ax2.YLim = [0.5 4.5]; Ax2.YLabel.String = 'Sleep State'; xlabel('Hrs'); title('Haar Approximation and Sleep State'); %% % The Haar approximation more clearly shows the relationship between the % sleep state and the heart rate data. You can assess this change by % looking at the correlation between the raw data and the sleep state time % series. corr(SS,HR) %% % Now compare the value of 0.56 with the correlation between the sleep % state data and the Haar approximation corr(SS,HaarHR) %% % The correlation has increased from 0.56 to 0.69. More advanced wavelet % analysis and modeling of this data is presented in Nason, von Sachs, & % Kroisandt (2000) and Nason, Sapatinas, & Sawczenko (2001). %% Digital Watermarking of Images % Watermarking is an important data protection tool. It is a passive % protection technique where a marker is covertly inserted in some data in % order to verify the authenticity or integrity of the data. Wavelet % techniques in general and the Haar transform in particular are frequently % employed in watermarking images. This example illustrates the use of the % Haar transform in watermarking an image and recovering the watermark. The % example employs a simple watermarking scheme chosen for ease of % illustration. In this scheme, the watermark is inserted in the % approximation coefficients at level 3. %% % Watermark an image of a Mandril with one of a honey badger. Read in the % Mandrill image. Resize it to 2048x2048 and display the result. coverIM = imread('mandrill.jpg'); coverIM = rgb2gray(coverIM); coverIM = imresize(im2double(coverIM),[2048 2048]); imagesc(coverIM); colormap gray; title('Original Image to Watermark'); axis off; axis square; %% % Obtain the Haar transform of the Mandrill image down to level 3. [LLorig,LHorig,HLorig,HHorig] = haart2(coverIM,3); imagesc(LLorig), title('Level 3 Haar Approximation'); axis off; axis square; %% % Read in watermark image and resize it. watermark = imread('honeybadger.jpg'); watermark = im2double(rgb2gray(watermark)); watermark = imresize(watermark,[2048 2048]); %% % Obtain the Haar transform of the watermark image down to level 3. [LLw,LHw,HLw,HHw] = haart2(watermark,3); imagesc(LLw); colormap gray; title('Level 3 Haar Approximation--Watermark'); axis off; axis square; %% % Add the honey badger watermark to the Mandril image by attenuating the % level-3 approximation coefficients of the watermark and inserting the % attenuated coefficients into the level-3 Mandrill approximation % coefficients. LLwatermarked = LLorig+0.0001*LLw; markedIM = ihaart2(LLwatermarked,LHorig,HLorig,HHorig); imagesc(markedIM); title('Watermarked Image') axis off; axis square; colormap gray; %% % The watermark (honey badger) is not visible in the watermarked image. % Because you know what algorithm was used to insert the watermark, you % can recover the watermark using the Haar transform. [LLr,LHr,HLr,HHr] = haart2(markedIM,3); LLmarked = (LLr-LLorig).*1e3; imagesc(LLmarked); title('Recovered Watermark'); colormap gray; axis off; axis square; %% References % Nason, G.P, von Sachs R. & Kroisandt, G. (2000) _Wavelet processes and % adaptive estimation of the evolutionary wavelet spectrum_. % J. R. Statist. Soc. Series B, 62, 271-292. % % Nason, G.P, Sapatinas, T. & Sawczenko, A. (2001) _Wavelet packet % modelling of infant sleep state using heart rate data_, Sankhya B, 63, % 199-217. % % Percival & Walden (2000) _Wavelet Methods for Time Series Analysis_, % Cambridge University Press.