c3ex10.m

来自「数字信号处理Matlab实现?中藕糯鞰atlab实现数字信号处理Matlab」· M 代码 · 共 28 行

%	Plots for Example 3-10
%
RC = input('Enter RC time constant of filter (square wave period = 1) ');
n_max = input('Enter vector of highest harmonic values desired (odd) ');
N = length(n_max);
t = 0:.005:1;
omega_0 = 2*pi;
for k = 1:N
	n = [-n_max(k):2:n_max(k)];
	L_n = length(n);
	nn = 2:L_n/2+1;
	sgn = (-1).^nn;		% This will make the 1,5,... terms + and the 3,7,...terms -
% 	Fourier coefficients for odd square wave; for even square wave, leave exp( ) off
	X_n = (2./(pi*abs(n))).*[fliplr(sgn) sgn].*exp(-j*0.5*pi*n);	
%	Input to filter
	x = X_n*exp(j*omega_0*n'*t);	% Rows of exponential matrix are versus time; 
					% columns are versus n; matrix multiply sums over n
%	Fourier coefficients of output
	Y_n = X_n.*(1./(1+j*n*omega_0*RC));
	y = Y_n*exp(j*omega_0*n'*t);	% The output
% 	Plot real part of y to get rid of small imaginary part due to computational error
	subplot(N,1,k),plot(t, real(y)), xlabel('t'), ylabel('approx. output'),...
	axis([0 1 -1.5 1.5]), text(.1,-.25, ['max. har. = ', num2str(n_max(k))])
	if k == 1
		title(['Filter output for RC = ', num2str(RC), ' seconds'])
	end
end