# 浅谈压缩感知（十四）：傅里叶矩阵与小波变换矩阵的MATLAB实现

1. 傅里叶矩阵及其MATLAB实现
2. 小波变换矩阵及其MATLAB实现

## 傅里叶矩阵及其MATLAB实现

dftmtx(N) is the N-by-N complex matrix of values around the unit-circle whose inner product with a column vector of length N yields the discrete Fourier transform of the vector. If X is a column vector of length N, then dftmtx(N)*X yields the same result as FFT(X); however, FFT(X) is more efficient.

The inverse discrete Fourier transform matrix is CONJ(dftmtx(N))/N.

% clc;clear;
N = 16;
X = randn(N,1);
dft_result1 = dftmtx(N)*X;
dft_result2 = fft(X);
% isEqual = all(dft_result1 == dft_result2);
err = norm(dft_result1(:)-dft_result2(:));
if err < 0.01
fprintf('dftmtx(N)*X yields the same result as FFT(X)');
else
fprintf('dftmtx(N)*X does not yield the same result as FFT(X)');
end

## 小波变换矩阵及其MATLAB实现

function [ ww ] = dwtmtx( N,wtype,wlev )
%DWTMTX Discrete wavelet transform matrix
%   This function generates the transform matrix ww according to input
%   parameters N,wtype,wlev .
%Detailed explanation goes here
%   N is the dimension of ww
%   wtype is the wavelet type
%   wlev is the number of decomposition level
%NOTE: The extension mode must be Periodization('per')
[h,g]= wfilters(wtype,'d');         %Decomposition low&high pass filter
L=length(h);                        %Filter length
h_1 = fliplr(h);                    %Flip matrix left to right
g_1 = fliplr(g);
loop_max = log2(N);
loop_min = double(int8(log2(L)))+1;
if wlev>loop_max-loop_min+1
fprintf('\nWaring: wlev is too big\n');
fprintf('The biggest wlev is %d\n',loop_max-loop_min+1);
wlev = loop_max-loop_min+1;
end
ww=1;
for loop = loop_max-wlev+1:loop_max
Nii = 2^loop;
p1_0 = [h_1 zeros(1,Nii-L)];
p2_0 = [g_1 zeros(1,Nii-L)];
p1 = zeros(Nii/2,Nii);
p2 = zeros(Nii/2,Nii);
for ii=1:Nii/2
p1(ii,:)=circshift(p1_0',2*(ii-1)+1-(L-1)+L/2-1)';
p2(ii,:)=circshift(p2_0',2*(ii-1)+1-(L-1)+L/2-1)';
end
w1=[p1;p2];
mm=2^loop_max-length(w1);
w=[w1,zeros(length(w1),mm);zeros(mm,length(w1)),eye(mm,mm)];
ww=ww*w;
clear p1;clear p2;
end

%The end!!!
end

%验证函数dwtmtx的正确性
clear all;close all;clc;
N = 2^7;
% 'db1' or 'haar', 'db2', ... ,'db10', ... , 'db45'
% 'coif1', ... , 'coif5'
% 'sym2', ... , 'sym8', ... ,'sym45'
% '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'
% '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'
wtype = 'rbio6.8';
wlev_max = wmaxlev(N,wtype);
if wlev_max == 0
fprintf('\nThe parameter N and wtype does not match!\n');
end
dwtmode('per');
for wlev = 1:wlev_max
ww = dwtmtx(N,wtype,wlev);
x = randn(1,N);
y1 = (ww*x')';
[y2,y2l] = wavedec(x,wlev,wtype);
y_err = sum((y1-y2).*(y1-y2));
fprintf('wlev = %d: y_err = %f\n',wlev,y_err);
end

posted @ 2015-12-25 15:51  AndyJee  阅读(10038)  评论(0编辑  收藏  举报