demofract.m

分形样条和小波的源程序,使用matlab代码

function demofract(alpha,beta,J)
% 	Generation of fractional-brownian-motion-like noise using a fractional
%	spline wavelet transform.
% 
% 	Author: Thierry Blu and Michael Unser, October 1999
% 	Biomedical Imaging Group, EPFL, Lausanne, Switzerland.
% 	This software is downloadable at http://bigwww.epfl.ch/
%
%	Example of parameters : J=12, alpha=0.33, beta=1.33
%
disp(' ')
disp('Generation of a fBm-like noise using orthogonal fractional splines')
disp(' ')

if nargin<3
	J=input(    'Choose J where  M=2^J is the length of input signal: ');
end
M=2^J;
if nargin<1
	alpha=input('Choose the degree of the spline (>-0.5): ');
end
if alpha<=-0.5
	disp(' ')
	disp(['The autocorrelation of the fractional splines exists only '...
	'for degrees strictly larger than -0.5!'])
	disp(' ')
	return
end
if nargin<2
	beta=input( 'Choose the power decrease of the inband power (typically, alpha+1): ');
end

disp(' ')

%
% Generate the wavelet transform: w=[wav1 wav2	... wavJ lowJ]
%

w=[];
s=M/2;
for j=1:J
	w=[w sqrt(12)*(rand(1,s)-0.5)*2^(j*(beta))];
	s=s/2;
end
w=[w sqrt(12)*(rand(1,2*s)-0.5)*2^(J*(beta))];

%
% Initializes wavelet filters and computes inverse wavelet transform
%

type='*ortho';
[FFTanalysisfilters,FFTsynthesisfilters]=FFTfractsplinefilters(M,alpha,type);
x=FFTwaveletsynthesis(w,FFTsynthesisfilters,J);

%
% Plotting of results
%

hold off
subplot(121)
plot(x)
xlabel('sample #')
title('Synthetized fBm-like noise')
a=axis;
axis([1 M a(3) a(4)])
subplot(122)
y=1/M*abs(fft(x)).^2;
y=[y((1+M/2):M) y(1:M/2)];
nu=abs([-M/2:(M/2-1)]/M);

semilogx(nu,10*log10(y))
a=axis;
axis([1/M 0.5 a(3) a(4)])
xlabel('normalized frequency')
ylabel('Fourier amplitude (dB)')
title('DSP of the synthetized noise')
set(gca,'XTick',[0.01 0.1 0.2 0.5],'XTickLabel',['0.01';' 0.1';' 0.2';' 0.5'])
hold on
h1=semilogx(nu,20*log10(nu.^(-beta)),'r');
h2=semilogx(nu,20*log10(sqrt(y(2+M/2))*(nu/nu(2+M/2)).^(-alpha-1)),'g');
legend([h1,h2],'inband power decrease','fractional spline decrease')
hold off

%
% Checks that the wavelet transform is reversible.
%

w0=FFTwaveletanalysis(x,FFTanalysisfilters,J);
x0=FFTwaveletsynthesis(w0,FFTsynthesisfilters,J);
disp(['Analysis-Synthesis SNR: ' ...
num2str(-round(20*log10(max(abs(x-x0))/max(abs(x))))) ' dB'])