c10ex13.m

数字信号处理Matlab实现?中藕糯鞰atlab实现数字信号处理Matlab实现

%	Study of weighting (window) functions - Example 10-13
%
N = input('Enter window size - points ');
NFFT = input('Enter number of FFT points ');
w = zeros(size(N));
K = menu('Weighting Function Options', 'rectangular',...
	'triangular',...
	'sine lobe',...
	'Hanning',...
	'sine cubed',...
	'sine to the fouth',...
	'Hamming',...
	'Blackman',...
	'Kaiser-Bessel',...
	'Gaussian');
zz = zeros(size(1:NFFT-N));
if K == 1
	w(1:N) = ones(1:N);
elseif K == 2
	for n = 1:N
		if n <= N/2+1
			w(n) = 2*(n-1)/N;
		else
			w(n) = 2*(N-n+1)/N;
		end
	end
elseif K == 3
	n = 1:N;
	w = sin(pi*(n-1)/N);
elseif K == 4
	n = 1:N;
	w = (sin(pi*(n-1)/N)).^2;
elseif K == 5
	n = 1:N;
	w = (sin(pi*(n-1)/N)).^3;
elseif K == 6
	n = 1:N;
	w = (sin(pi*(n-1)/N)).^4;
elseif K == 7
	n = 1:N;
	w = 0.54 - 0.46*cos(2*pi*(n-1)/N);
elseif K == 8
	n = 1:N;
	w = 0.42 - 0.50*cos(2*pi*(n-1)/N) + 0.08*cos(4*pi*(n-1)/N);
elseif K == 9
	alpha = input('Enter value for parameter alpha ');
	arg = zeros(size(N));
	for n = 1:N
		arg(n) = pi*alpha*sqrt(1-(2*(n-N/2-1)/N)^2);
	end
	w = besselI(0,arg)/besselI(0,pi*alpha);
elseif K == 10
	alpha = input('Enter value for parameter alpha ');
	n = 1:N;
	w = exp(-.5*(2*alpha*(n-N/2-1)/N).^2);
end
w = [w zz];
W = fft(w);
subplot(2,1,1),stem(w), axis([1 2*N 0 2]),xlabel('n'),ylabel('w(n)'),...
	if K == 1
	text(15, 1.5, ['Rect. window; L = ',num2str(N)])
	elseif K == 2
	text(15, 1.5, ['Triang. window; L = ',num2str(N)])
	elseif K == 3
	text(15, 1.5, ['Sine window; L = ',num2str(N)])
	elseif K == 4
	text(15, 1.5, ['Hanning window; L = ',num2str(N)])
	elseif K == 5
	text(15, 1.5, ['Sine^3 window; L = ',num2str(N)])
	elseif K == 6
	text(15, 1.5, ['Sine^4 window; L = ',num2str(N)])
	elseif K == 7
	text(15, 1.5, ['Hamming window; L = ',num2str(N)])
	elseif K == 8
	text(15, 1.5, ['Blackman window; L = ',num2str(N)])
	elseif K == 9
	text(15, 1.5, ['Kaiser-Bessel window; L = ',num2str(N)])
	elseif K == 10
	text(15, 1.5, ['Gaussian window; L = ',num2str(N)])
	end
subplot(2,1,2),plot(20*log10((abs(W)+.00000001)/max(abs(W)))),...
	 axis([1 NFFT -80 0]),xlabel('n'),ylabel('FFT of window - dB')