%% Exercises from Linear algebra, signal processing, and wavelets. A unified approach.\\ Matlab version
% 
% Author: \O yvind Ryan
% 
% Date: Feb 12, 2017
% 
% % Externaldocuments: applinalg
% % Mapping from exercise labels to numbers: label2numbers = {'playdwtex2': '5.28', 'playdwtex3': '5.35', 'fdwtex3': '5.36', 'fdwtex2': '5.29', 'fdwtex1': '5.11', 'oblig2ex1': '5.39', 'simpledwtexample': '5.9', 'ex050316': '5.25', 'example:cascade_pwc_pwl': '5.43', 'listendifference': '5.42', 'playdwtex': '5.10'}
% 
% 
% Table of contents:
% 
%  Sound and Fourier series 
%  Digital sound and Discrete Fourier analysis 
%  Operations on digital sound: digital filters 
%  Symmetric filters and the DCT 
%  Motivation for wavelets and some simple examples 
%      Example 5.9: Computing the DWT by hand 
%      Example 5.10: DWT on sound 
%      Example 5.11: DWT on the samples of a mathematical function 
%      Exercise 5.25: Computing the DWT of a simple vector 
%      Example 5.28: DWT on sound 
%      Example 5.29: DWT on the samples of a mathematical function 
%      Example 5.35: DWT on sound 
%      Example 5.36: DWT on the samples of a mathematical function 
%      Exercise 5.39: Implement finding $\psi$ with vanishing moments 
%      Exercise 5.42: Listening experiments 
%      Example 5.43: Implementing the cascade algorithm 
%  The filter representation of wavelets 
%  Constructing interesting wavelets 
%  The polyphase representation and wavelets 
%  Digital images 
%  Using tensor products to apply wavelets to images 
%  Appendix: Basic Linear Algebra 
%  Appendix: Signal processing and linear algebra: a translation guide 
% 
% 
% 
%% Sound and Fourier series
%% Digital sound and Discrete Fourier analysis
%% Operations on digital sound: digital filters
%% Symmetric filters and the DCT
%% Motivation for wavelets and some simple examples
% 
% 
% % --- begin exercise ---
% 
%% Example 5.9: Computing the DWT by hand
% 
% % keywords = ipynbwavexamples; mwavexamples
% 
% In some cases, the DWT can be computed by hand, keeping in mind its definition as a change of coordinates. 
% As an example, consider the simple vector $\bm{x}$ of length $2^{10}=1024$ defined by
% 
% $$\begin{equation*} x_n=\begin{cases} 1 & \text{for } n< 512 \\ 0 & \text{for } n\geq 512,\end{cases} \end{equation*}$$
% and let us compute the $10$-level DWT of this vector by first visualizing the function with these coordinates. 
% Since $m=10$ here, we should view $\bm{x}$ as coordinates in the basis $\vphi_{10}$ of a function $f(t)\in V_{10}$. 
% This is $f(t)=\sum_{n=0}^{511} \phi_{10,n}$, and since $\phi_{10,n}$ is supported on
% $[2^{-10}n,2^{-10}(n+1))$, the support of $f$ has width $512\times 2^{-10}=1/2$ (512 translates, each with width $2^{-10}$). 
% Moreover, since $\phi_{10,n}$ is $2^{10/2}=2^5=32$ on $[2^{-10}n,2^{-10}(n+1))$ and $0$ elsewhere, it is clear that
% 
% $$\begin{equation*} f(t) = \begin{cases} 32 & \text{for } 0\leq t<1/2 \\ 0 & \text{for } 0 t\geq 1/2. \end{cases} \end{equation*}$$
% This is by definition a function in $V_1$: $f$ must in fact be a multiplum of $\phi_{1,0}$, since this also is supported on $[0,1/2)$. 
% We can thus write $f(t)=c\phi_{1,0}(t)$ for some $c$. We can find $c$ by setting $t=0$. This gives that
% $32=2^{1/2}c$ (since $f(0)=32$, $\phi_{1,0}(0)=2^{1/2}$), so that $c=32/\sqrt{2}$. 
% This means that $f(t)=\frac{32}{\sqrt{2}}\phi_{1,0}(t)$, $f$ is in $V_1$, and with coordinates $(32/\sqrt{2},0,\ldots,0)$ in $\vphi_1$.
% 
% When we run a $10$-level DWT we make a change of coordinates from $\vphi_{10}$ to $(\vphi_0,\vpsi_0,\cdots,\vpsi_{9})$. 
% The first $9$ levels give us the coordinates in $(\vphi_1,\vpsi_1,\vpsi_2,\ldots,\vpsi_9)$, and these are $(32/\sqrt{2},0,\ldots,0)$ from what we showed.
% It remains thus only to perform the last level in the DWT, i.e. perform the change of coordinates from $\vphi_1$ to $(\vphi_0,\vpsi_0)$. 
% Since $\phi_{1,0}=\frac{1}{\sqrt{2}}(\phi_{0,0}+\psi_{0,0})$, so that we get
% 
% $$\begin{equation*} f(t)=\frac{32}{\sqrt{2}}\phi_{1,0}(t)=\frac{32}{\sqrt{2}}\frac{1}{\sqrt{2}}(\phi_{0,0}+\psi_{0,0})=16\phi_{0,0}+16\psi_{0,0}. \end{equation*}$$
% From this we see that the coordinate vector of $f$ in $(\vphi_0,\vpsi_0,\cdots,\vpsi_{9})$, i.e. the 10-level DWT of $\bm{x}$,  is 
% $(16,16,0,0,\ldots,0)$. Note that here $V_0$ and $W_0$ are both $1$-dimensional, since $V_{10}$ was assumed to
% be of dimension $2^{10}$ (in particular, $N=1$).
% 
% It is straightforward to verify what we found using the algorithm above:
% 
DWTImpl([ones(512,1); zeros(512,1)], 10, 'Haar')
% The reason why the method from this example worked was that the vector we started with had a simple representation in the wavelet basis, actually it equaled
% the coordinates of a basis function in $\vphi_1$. Usually this is not the case, and our only possibility then is to run the DWT on a computer.
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 5.10: DWT on sound
% 
% % keywords = ipynbwavexamples; mwavexamples
% 
% Let us plot the samples of our audio sample file, and compare them with the first order DWT. 
% 
[x,fs]=audioread('sounds/castanets.wav');
nvals = (0:2^17)';
figure();
plot(nvals(1:100:end),x(1:100:2^17,1),'k-')
% 
newx=x(1:2^17,1);
newx = DWTImpl(newx,1,'Haar');
figure();
plot(nvals(1:100:end),newx(1:100:end),'k-')
% The first part of the DWT plot represents the low resolution part, the second the detail. 
% 
% Since $\phi(2^mt-n)\in V_m$ oscillates more quickly than $\phi(t-n)\in V_0$,
% one is lead to believe that coefficients from lower order resolution spaces correspond to lower frequencies. 
% The functions $\phi_{m,n}$ do not correspond to pure tones in the setting of wavelets, however,  
% but let us nevertheless listen to sound from the different resolution spaces. 
% The code base includes a function |forw_comp_rev_DWT| which runs an $m$-level DWT on the first samples 
% of the audio sample file, extracts the detail or the low-resolution approximation, and runs an IDWT to reconstruct the sound. 
% Since the returned values may lie outside the legal range $[-1,1]$, the values are normalized at the end.
% 
% 
% To listen to the low-resolution approximation, write
% 
m = 1;
% 
[x, fs] = forw_comp_rev_DWT(m, 'Haar');
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% It is instructive to run this code for different values of $m$. 
% For $m=2$ we clearly hear a degradation in the sound. 
% For $m=4$ and above most of the sound is unrecognizable, since too much of the detail is omitted. 
% To be more precise, when listening to the sound by throwing away detail from $W_0$, $W_1$,...,$W_{m-1}$, we are left with a $2^{-m}$ share of the data. 
% 
% Let us also consider the detail. For $m=1$ this can be played as follows
% 
[x, fs] = forw_comp_rev_DWT(1, 'Haar', 0);
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% and plotted as follows
% 
plot((1:100:2^17)', x(1:100:end),'k-')
% We see that the detail is quite significant, so that the first order wavelet approximation does not give a very good approximation. 
% For $m=2$ the detail can be played as follows
% 
[x, fs] = forw_comp_rev_DWT(2, 'Haar', 0);
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% and plotted as follows
% 
plot((1:100:2^17)', x(1:100:end),'k-')
% The error is larger when two levels of the DWT are performed, as one would suspect. 
% It is also seen that the error is larger in the part of the file where there are bigger variations. 
% Since more and more information is contained in the detail components as we increase $m$, we here see the opposite effect: 
% The sound gradually improves in quality as we increase $m$. 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% 
% The previous example illustrates that wavelets as well may be used to perform operations on sound. As we will see later, however, our main application for wavelets
% will be images, where they have found a more important role than for sound. Images typically display variations which are less abrupt than the ones found in sound.
% Just as the functions above had smaller errors in the corresponding resolution spaces than the sound had, images are thus more suited for
% for use with wavelets. The main idea behind why wavelets are so useful comes from the fact that the detail, i.e.,
% wavelet coefficients corresponding to the spaces $W_k$, are often very small. After a DWT one is therefore often left with a couple of significant coefficients,
% while most of the coefficients are small. The approximation from $V_0$ can be viewed as a good approximation, even though it contains much less
% information. This gives another reason why wavelets are popular for images: Detailed images can be very large, but when they are downloaded to a
% web browser, the browser can very early show a low-resolution of the image, while waiting for the rest of the details in the image to be downloaded. When we later
% look at how wavelets are applied to images, we will need to handle one final hurdle, namely that images are two-dimensional.
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 5.11: DWT on the samples of a mathematical function
% 
% % keywords = ipynbwavexamples; mwavexamples
% 
% Above we plotted the DWT coefficients of a sound, as well as the detail/error. We can also experiment with samples generated from a mathematical function.
% The following code computes input vectors corresponding to a piecewise constant, a trigonometric, and a piecewise linear function. 
% 
N = 2^10;
nvals = (linspace(0,N-1,N))';
fpwc = ones(N,1);
fpwc(1:(N/2)) = fpwc(1:(N/2))*0.75;
fpwc((N/2+1):end) = fpwc((N/2+1):end)*0.25;
ftrig = 0.5+0.5*cos(2*pi*nvals/N);
flin =  nvals*2/(N-1);
flin = 1-abs(1-flin);
% Finally, the following code compares the computed details with the entries in the vectors. 
% 
[detail1, detail2] = findetailm6m8(fpwc, 'Haar');
plot(nvals, abs(fpwc),'r.', nvals, abs(detail1), 'g.',nvals,abs(detail2),'b.')
% 
[detail1, detail2] = findetailm6m8(ftrig, 'Haar');
plot(nvals,abs(ftrig),'r.',nvals,abs(detail1),'g.',nvals,abs(detail2),'b.')
% 
[detail1, detail2] = findetailm6m8(flin, 'Haar');  
plot(nvals,abs(flin),'r.',nvals,abs(detail1),'g.',nvals,abs(detail2),'b.')
% In these cases, we see that we require large $m$ before the detail/error becomes significant. We see also that there is no error for the square wave. The reason is
% that the square wave is a piecewise constant function, so that it can be represented exactly by the $\phi$-functions. For the other functions, however, this is not
% the case, so we here get an error.
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Exercise 5.25: Computing the DWT of a simple vector
% 
% % keywords = ipynbwavexamples; mwavexamples
% 
% Suppose that we have the vector $\bm{x}$ with length $2^{10}=1024$, defined by $x_n=1$ for $n$ even, $x_n=-1$ for $n$ odd. 
% What will be the result if you run a 10-level DWT on $\bm{x}$?
% Use the function |DWTImpl| to verify what you have found.
% 
% % --- begin hint in exercise ---
% 
% *Hint.* We defined $\psi$ by $\psi(t)=(\phi_{1,0}(t)-\phi_{1,1}(t))/\sqrt{2}$. From this connection it follows that 
% $\psi_{9,n}=(\phi_{10,2n}-\phi_{10,2n+1})/\sqrt{2}$, and thus $\phi_{10,2n}-\phi_{10,2n+1}=\sqrt{2}\psi_{9,n}$. 
% Try to couple this identity with the alternating sign you see in $\bm{x}$.
% 
% % --- end hint in exercise ---
% 
% 
% % --- begin solution of exercise ---
% *Solution.* The vector $\bm{x}$ is the coordinate vector of the function $f(t)=\sum_{n=0}^{1023} (-1)^n \phi_{10,n}$ in the basis $\vphi_{10}$ for $V_{10}$. 
% Since $\phi_{10,2n}-\phi_{10,2n+1}=\sqrt{2}\psi_{9,n}$, we can write $f(t)=\sum_{n=0}^{1023} \sqrt{2}\psi_{9,n}$. 
% Since a 10-level-DWT gives as a result the coordinate vector of $f$ in
% 
% $$\begin{equation*} (\vphi_0,\vpsi_0,\vpsi_1,\vpsi_2,\vpsi_3,\vpsi_4,\vpsi_5,\vpsi_6,\vpsi_7,\vpsi_8,\vpsi_9), \end{equation*}$$
% (the DWT is nothing but the change of coordinates from $\vphi_{10}$ to this basis), and since $f(t)=\sum_{n=0}^{1023} \sqrt{2}\psi_{9,n}$, it is clear that the
% coordinate vector of $f$ in this basis has $\sqrt{2}$ in the second part (the $\vpsi_9$-coordinates), and $0$ elsewhere.
% The 10-level DWT of $\bm{x}$ therefore gives the vector of length 1024 which is $0$ on the first half, and equal to $\sqrt{2}$ on the second half. 
% $m=10$ is here arbitrarily chosen: The result would have been the same for $m=1,m=2,$ and so on.
% The following code verifies the result:
% 
DWTImpl(repmat([1; -1],512,1), 10, 'Haar')
% 
% % --- end solution of exercise ---
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 5.28: DWT on sound
% 
% % keywords = ipynbwavexamples; mwavexamples
% 
% Using the new kernels, let us plot the listen to the new low resolution approximations, 
% as well as plot and listen to the detail, as we did in Example 5.10. 
% First we listen to the low-resolution approximation.
% 
% 
m = 1;
% 
[x, fs] = forw_comp_rev_DWT(m, 'pwl0');
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% There is a new and undesired effect when we increase $m$ here: The castanet sound seems to grow strange. 
% The sounds from the castanets are perhaps the sound with the highest frequencies. 
% 
% Now for the detail. For $m=1$ this can be played as follows
% 
[x, fs] = forw_comp_rev_DWT(1, 'pwl0', 0);
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% and plotted as follows
% 
plot((1:100:2^17)', x(1:100:end),'k-')
% For $m=2$ the detail can be played as follows
% 
[x, fs] = forw_comp_rev_DWT(2, 'pwl0', 0);
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% and plotted as follows
% 
plot((1:100:2^17)', x(1:100:end),'k-')
% When comparing with Example 5.10 we see much of the same, but it seems here that the error is bigger than before. 
% In the next section we will try to explain why this is the case, and construct another wavelet based on piecewise linear functions which remedies this.
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 5.29: DWT on the samples of a mathematical function
% 
% % keywords = ipynbwavexamples; mwavexamples
% 
% Let us also repeat Example 5.11, where we plotted the detail/error at different resolutions for the samples of a mathematical function. 
% 
% We first compute new input vectors corresponding to a piecewise constant, a trigonometric, and a piecewise linear function. 
% 
N = 2^10 + 1;
nvals = (linspace(0,N-1,N))';
fpwl = ones(N,1);
fpwl(1:(N+1)/2) = fpwl(1:(N+1)/2)*0.75;
fpwl(((N+1)/2+1):end) = fpwl(((N+1)/2+1):end)*0.25;
ftrig = 0.5+0.5*cos(2*pi*nvals/N);
flin =  nvals*2/(N-1);
flin = 1-abs(1-flin);
% The following code compares the computed details with the entries in the new vectors. 
% 
[detail1, detail2] = findetailm6m8(fpwl, 'pwl0');
plot(nvals,abs(fpwl),'r.',nvals,abs(detail1),'g.',nvals,abs(detail2),'b.')
% 
[detail1, detail2] = findetailm6m8(ftrig, 'pwl0');
plot(nvals,abs(ftrig),'r.',nvals,abs(detail1),'g.',nvals,abs(detail2),'b.')
% 
f =  nvals*2/(N-1); 
f = 1-abs(1-f);
[detail1, detail2] = findetailm6m8(flin, 'pwl0');
plot(nvals,abs(flin),'r.',nvals, abs(detail1),'g.', nvals, abs(detail2),'b.')
% With the square wave we see now that there is an error. The reason is that a piecewise constant function can not be represented exactly by piecewise linear
% functions, due to discontinuity. For the second function we see that there is no error. The reason is that this function is piecewise linear, so there is no
% error when we represent the function from the space $V_0$. 
% With the third function, however, we see an error.
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 5.35: DWT on sound
% 
% % keywords = ipynbwavexamples; mwavexamples
% 
% Using the new kernels, let us also here listen to the low resolution approximations and the detail.  
% First the low-resolution approximation:
% 
% 
m = 1;
% 
[x, fs] = forw_comp_rev_DWT(m, 'pwl2');
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% The new, undesired effect in the castanets from Example 5.28 now seem to be gone. 
% The detail for $m=1$ this can be played as follows
% 
[x, fs] = forw_comp_rev_DWT(1, 'pwl2', 0);
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% and plotted as follows
% 
vals=(1:100:2^17)';
plot(vals, x(1:100:end),'k-')
% For $m=2$ the detail can be played as follows
% 
[x, fs] = forw_comp_rev_DWT(2, 'pwl2', 0);
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% and plotted as follows
% 
vals=(1:100:2^17)';
plot(vals, x(1:100:end),'k-')
% Again, when comparing with Example 5.10 we see much of
% the same. It is difficult to see an improvement from this figure. However, this figure also clearly shows a smaller error than the piecewise linear wavelet. 
% A partial explanation is that the wavelet we now have constructed has two vanishing moments, while the other had not.
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 5.36: DWT on the samples of a mathematical function
% 
% % keywords = ipynbwavexamples; mwavexamples
% 
% Let us also repeat Exercise 5.11 for our alternative wavelet, where we plotted the detail/error at different resolutions, for the samples of a
% mathematical function. 
% The following code compares the computed details with the entries in the vectors. 
% 
[detail1, detail2] = findetailm6m8(fpwl, 'pwl2');
plot(nvals,abs(fpwl),'r.', nvals,abs(detail1),'g.', nvals,abs(detail2),'b.')
% 
[detail1, detail2] = findetailm6m8(ftrig, 'pwl2');
plot(nvals,abs(ftrig),'r.',nvals,abs(detail1),'g.',nvals,abs(detail2),'b.')
% 
[detail1, detail2] = findetailm6m8(flin, 'pwl2');
plot(nvals,abs(flin),'r.',nvals, abs(detail1),'g.',nvals, abs(detail2),'b.')
% Again for the square wave there is an error, which seems to be slightly lower than for the previous wavelet. For the second function we see that there is no error,
% as before. The reason is the same as before, since the function is piecewise linear. With the third function there is an error. The error seems to be slightly
% lower than for the previous wavelet, which fits well with the fact that this new wavelet has a bigger number of vanishing moments.
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Exercise 5.39: Implement finding $\psi$ with vanishing moments
% 
% % keywords = ipynbwavexamples; mwavexamples
% 
% In the previous exercise we ended up with a lot of calculations to find $\alpha,\beta$ in 
% Equation (5.38). Let us try to make a program which does
% this for us, and which also makes us able to generalize the result.
% 
% 
% *a)* Define
% 
% $$\begin{align*} a_k &= \int_{-1}^1 t^k(1-|t|)dt, & b_k &= \int_0^2 t^k(1-|t-1|)dt, & e_k &= \int_0^1 t^k(1-2|t-1/2|)dt, \end{align*}$$
% for $k\geq 0$. Explain why finding $\alpha,\beta$ so that we have two vanishing moments in Equation 
% (5.38) 
% is equivalent to solving the following equation:
% 
% $$\begin{equation*} \begin{pmatrix} a_0 & b_0 \\  a_1 & b_1 \end{pmatrix} \begin{pmatrix} \alpha \\ \beta \end{pmatrix} = \begin{pmatrix} e_0 \\ e_1 \end{pmatrix} \end{equation*}$$
% Write a program which sets up and solves this system of equations, and use this program to verify the values for $\alpha,\beta$ we previously have found.
% 
% % --- begin hint in exercise ---
% 
% *Hint.* you can integrate functions in Matlab with the function |quad|
% . As an example, the function $\phi(t)$, which is nonzero only on $[-1,1]$, can be integrated as follows:
% 
%  quad(@(t)t.^k.*(1-abs(t)),-1,1)
% 
% % --- end hint in exercise ---
% 
% 
% % --- begin solution of exercise ---
% *Solution.* In order for $\psi$ to have vanishing moments we must have that $\int\hat\psi(t)dt=\int t\hat\psi(t)dt=0$
% Substituting $\hat\psi=\psi-\alpha \phi_{0,0}-\beta \phi_{0,1}$ we see that, for $k=0,1$,
% 
% $$\begin{equation*} \int t^k\left(\alpha \phi_{0,0}+\beta \phi_{0,1}\right)dt = \int t^k\psi(t)dt. \end{equation*}$$
% The left hand side can here be written
% 
% $$\begin{align*} \int t^k\left(\alpha \phi_{0,0}+\beta \phi_{0,1}\right)dt &= \alpha\int t^k\phi_{0,0}dt + \beta\int t^k\phi_{0,1}(t)dt \\  &=\alpha\int_{-1}^1 t^k(1-|t|)dt + \beta\int_0^2 t^k(1-|t-1|)dt =\alpha a_k + \beta b_k. \end{align*}$$
% The right hand side is
% 
% $$\begin{equation*}\int t^k\psi(t)dt = \int t^k\phi_{1,1}(t)dt = \int_0^1 (1-2|t-1/2|)dt=e_k.\end{equation*}$$
% The following program sets up the corresponding equation systems, and solves it by finding $\alpha,\beta$. 
% 
A = zeros(2);
b = zeros(2, 1);
for k = 0:1
A(k + 1, :) = [quad(@(t)t.^k.*(1 - abs(t)), -1, 1)...
quad(@(t)t.^k.*(1 - abs(t - 1)), 0, 2)];
b(k + 1) = quad(@(t)t.^k.*(1 - 2*abs(t - 1/2)), 0, 1);
end
A\b;
% 
% % --- end solution of exercise ---
% 
% *b)* The procedure where we set up a matrix equation in a) allows for generalization to more vanishing moments. Define
% 
% $$\hat\psi=\psi_{0,0}-\alpha \phi_{0,0}-\beta \phi_{0,1} -\gamma\phi_{0,-1}-\delta\phi_{0,2}.$$
% We would like to choose $\alpha,\beta,\gamma,\delta$ so that we have $4$ vanishing moments. Define also
% 
% $$\begin{align*} g_k &= \int_{-2}^0 t^k(1-|t+1|)dt, & d_k &= \int_1^3 t^k(1-|t-2|)dt \end{align*}$$
% for $k\geq 0$. Show that $\alpha,\beta,\gamma,\delta$ must solve the equation
% 
% $$\begin{equation*} \begin{pmatrix} a_0 & b_0 & g_0 & d_0 \\  a_1 & b_1 & g_1 & d_1 \\  a_2 & b_2 & g_2 & d_2 \\  a_3 & b_3 & g_3 & d_3 \end{pmatrix} \begin{pmatrix} \alpha \\ \beta \\ \gamma\\ \delta \end{pmatrix} = \begin{pmatrix} e_0 \\ e_1 \\ e_2 \\ e_3 \end{pmatrix}, \end{equation*}$$
% and solve this with your computer.
% 
% 
% % --- begin solution of exercise ---
% *Solution.* Similarly to a), Equation (5.47) gives that
% 
% $$\begin{equation*} \int t^k\left(\alpha \phi_{0,0}+\beta \phi_{0,1} +\gamma\phi_{0,-1}+\delta\phi_{0,2}\right)dt = \int t^k\psi(t)dt.  \end{equation*}$$
% The corresponding equation system is deduced exactly as in a). 
% The following program sets up the corresponding equation systems, and solves it by finding $\alpha,\beta,\gamma,\delta$.
% 
A = zeros(4);
b = zeros(4, 1);
for k = 0:3
A(k + 1, :) = [quad(@(t)t.^k.*(1 - abs(t)), - 1, 1)...
quad(@(t)t.^k.*(1 - abs(t - 1)), 0, 2)...
quad(@(t)t.^k.*(1 - abs(t + 1)), -2, 0)...
quad(@(t)t.^k.*(1 - abs(t - 2)), 1, 3)];
b(k + 1) = quad(@(t)t.^k.*(1 - 2*abs(t - 1/2)), 0, 1);
end
coeffs=A\b;
% 
% % --- end solution of exercise ---
% 
% *c)* Plot the function defined by (5.47), which you found in b).
% 
% % --- begin hint in exercise ---
% 
% *Hint.* If |t| is the vector of $t$-values, and you write
% 
%  (t >= 0).*(t <= 1).*(1-2*abs(t-0.5))
% you get the points $\phi_{1,1}(t)$.
% 
% % --- end hint in exercise ---
% 
% 
% % --- begin solution of exercise ---
% *Solution.* The function $\hat\psi$ now is supported on $[-2,3]$, and can be plotted as follows:
% 
t=linspace(-2,3,100);
plot(t, ( t>= 0).*(t <= 1).*(1-2*abs(t - 0.5)) ...
-coeffs(1)*(t >= -1).*(t <= 1).*(1 - abs(t))...
-coeffs(2)*(t >= 0).*(t <= 2).*(1 - abs(t - 1))...
-coeffs(3)*(t >= -2).*(t <= 0).*(1 - abs(t + 1))...
-coeffs(4)*(t >= 1).*(t <= 3).*(1 - abs(t - 2)))
% 
% % --- end solution of exercise ---
% 
% *d)* Explain why the coordinate vector of $\hat{\psi}$ in the basis $(\vphi_0,\vpsi_0)$ is
% 
% $$\begin{equation*}[\hat{\psi}]_{(\vphi_0,\vpsi_0)}=(-\alpha,-\beta,-\delta,0,\ldots,0-\gamma)\oplus(1,0,\ldots,0).\end{equation*}$$
% 
% % --- begin hint in exercise ---
% 
% *Hint.* The placement of $-\gamma$ may seem a bit strange here, and has to with that $\phi_{0,-1}$
% is not one of the basis functions $\{\phi_{0,n}\}_{n=0}^{N-1}$. However, we have that $\phi_{0,-1}=\phi_{0,N-1}$, i.e. $\phi(t+1)=\phi(t-N+1)$, since we always
% assume that the functions we work with have period $N$.
% 
% % --- end hint in exercise ---
% 
% *e)* Sketch a more general procedure than the one you found in b)., which can be used to find wavelet bases where we have even more vanishing moments.
% 
% 
% % --- begin solution of exercise ---
% *Solution.* If we define
% 
% $$\begin{equation*} \hat\psi=\psi_{0,0} - \sum_{k=0}^K \left( \alpha_k \phi_{0,-k} -\beta_k\phi_{0,k+1}\right), \end{equation*}$$
% we have $2k$ unknowns. These can be determined if we require $2k$ vanishing moments.
% 
% % --- end solution of exercise ---
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Exercise 5.42: Listening experiments
% 
% % keywords = ipynbwavexamples; mwavexamples
% 
% Run the function |forw_comp_rev_DWT| for different $m$ for the Haar wavelet, 
% the piecewise linear wavelet, and the alternative piecewise linear wavelet, but listen to the detail components 
% $W_0\oplus W_1\oplus\cdots\oplus W_{m-1}$ instead. 
% Describe the sounds you hear for different $m$, and try to explain why the sound seems to get louder when you increase $m$.
% 
% 
% % --- begin solution of exercise ---
% *Solution.* The following code can be used:
% 
[x, fs] = forw_comp_rev_DWT(m, 'Haar');
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% 
[x, fs] = forw_comp_rev_DWT(m, 'pwl0');
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% 
[x, fs] = forw_comp_rev_DWT(m, 'pwl2');
playerobj = audioplayer(x, fs);
playblocking(playerobj)
% 
% % --- end solution of exercise ---
% 
% % --- end exercise ---
% 
% 
% 
% 
% % --- begin exercise ---
% 
%% Example 5.43: Implementing the cascade algorithm
% 
% % keywords = ipynbwavexamples; mwavexamples
% 
% One thing should be noted in the function |cascade_alg|. 
% As the scaling function of the piecewise linear wavelet, it may be that the function is nonzero for small, negative values. 
% If we plot the function over $[0,N]$, we would see two disconnected segments - one to the left, and one to the right. 
% In the code we shift the values so that the graph appears as one connected segment. 
% 
% We will use $a=-2$ and $b=6$ in what follows, since $[-2,6]$ will turn out to contain all supports. We will also use $m=10$ levels in the cascade algorithm. 
% The following code then runs the cascade algorithm for the three wavelets we have considered, to reproduce all previous scaling functions and mother wavelets. 
% 
cascade_alg(10, -2, 6, 'Haar', 1, 0)
cascade_alg(10, -2, 6, 'Haar', 0, 0)
% 
cascade_alg(10, -2, 6, 'pwl0', 1, 0)
cascade_alg(10, -2, 6, 'pwl0', 0, 0)
% 
cascade_alg(10, -2, 6, 'pwl2', 1, 0)
cascade_alg(10, -2, 6, 'pwl2', 0, 0)
% In the above code you are encouraged to also try the other wavelets we encounter later in these notes. You may then also have to change the support interval $[a,b]$ accordingly.
% 
function [detail1, detail2]=findetailm6m8(x, wave_name)
N = size(x,1);
m6zeroout = N;
for mres = 0:5
m6zeroout = ceil(m6zeroout/2);
end
m8zeroout = m6zeroout;
m8zeroout = ceil(m8zeroout/2);
m8zeroout = ceil(m8zeroout/2);
m=6;
detail1=x;
detail1 = DWTImpl(detail1, m, wave_name);
detail1(1:m6zeroout) = 0;
detail1 = IDWTImpl(detail1, m, wave_name);
m=8;
detail2=x;
detail2 = DWTImpl(detail2, m, wave_name);
detail2(1:m8zeroout) = 0;
detail2 = IDWTImpl(detail2, m, wave_name);
end
% 
% % --- end exercise ---
% 
% 
%% The filter representation of wavelets
%% Constructing interesting wavelets
%% The polyphase representation and wavelets
%% Digital images
%% Using tensor products to apply wavelets to images
%% Appendix: Basic Linear Algebra
%% Appendix: Signal processing and linear algebra: a translation guide

