Path: news.mathworks.com!not-for-mail
From: <HIDDEN>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Negative indices in arrays
Date: Sun, 14 Jun 2009 02:38:01 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 220
Message-ID: <h11nq9$73l$1@fred.mathworks.com>
References: <h11ecd$mps$1@fred.mathworks.com> <h11jj9$fd9$1@fred.mathworks.com>
NNTP-Posting-Host: webapp-05-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1244947081 7285 172.30.248.35 (14 Jun 2009 02:38:01 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sun, 14 Jun 2009 02:38:01 +0000 (UTC)
Xref: news.mathworks.com comp.soft-sys.matlab:547266

Alright, welll. For an assignment I'm trying to calculate and plot the values for n for when they are zero and negative for the magnitude and phase plot at the very end of the Prob2.m for Dn. The negative indices are giving me the most difficulties :

Prob2.m:

clc
clear all
close all
format compact
%format short eng

% Create a waveform from t=Tbegin to t=Tend
Tbegin = 0.0;
Tend = 5;

N = 2^15; Nm1=N-1; Np1=N+1; Nd2=N/2;

t=[0:Nm1]';
t=Tbegin + t*Tend/Nm1;
tsample = t(2)-t(1);

T0 = (Tend - Tbegin);
f0 = 1/T0;
w0 = 2*pi*f0;
fnyq = f0/2; wnyq=2*pi*fnyq;

NNN=100;

t1 = 2;
x = ((-2*t).*(rect_fcn((t-1.5)/3)))+(rect_fcn((t-4)/2)) ;
%x(t<=t1) = t(t<=t1);

% evaluate the analytic expression
n = 1:NNN;
Dn = ((((exp(-j*n*(6/5)*pi)).*(2-(2*j*n*(6/5)*pi)))-2)./((n.^2)*(4/5)*(pi^2)))+(((exp(-j*n*2*pi))-(exp(-j*n*(6/5)*pi)))./(-j*n*2*pi));

[a0,a,b]=fourier_series_coef(t,x,NNN);
% Uncomment the following line to see the coeficients values.
%a0, [a b];

Dn2a = 2*real(Dn);
Dn2b = -2*imag(Dn);
% had to make N=2^10 larger at the beginning to get better
% numerical integration results
NNNd = 10;
nnd=(1:NNNd);

disp(' a_n coef ')
[nnd' a(1:NNNd) Dn2a(1:NNNd)']
disp(' ')
disp(' b_n coef ')
[nnd' b(1:NNNd) Dn2b(1:NNNd)']
disp(' ')
disp(' Dn')
[nnd' Dn(1:NNNd).']
disp(' ')
disp(' |Dn|')
[nnd' abs(Dn(1:NNNd))']
disp(' ')
disp(' arg(Dn)')
[nnd' angle(Dn(1:NNNd))']

% use the FFT to compute ~Dn
XDn = fft(x)/N;
disp(' ')
disp(' ')
disp(' analytic Dn FFT Dn')
[Dn(1:NNNd).' XDn(2:NNNd+1)]

% Generate a waveform using the Fourier Series
xg = fourier_series_gen(t,T0,a0,a,b);
subplot(2,1,1),plot(t,x), grid on
ylabel('\it{x(t)}')
title('One Period of \it{x(t)}');

subplot(2,1,2), plot(t,xg), grid on
xlabel('t (sec)')
ylabel('\it{x^''(t)}');
titlestr = ['Fourier Series Approximation to {\itx(t)} (' ...
num2str(NNN) ' terms)'];
title(titlestr);

% Compare using partial sums of 5, 15 and 50 coeficients

Np=3;
tp=[0:(Np*N-1)]';
tp=Tbegin + tp*Np*Tend/(Np*N-1);

% 5
Ncoef = 5;
[a0,a,b]=fourier_series_coef(t,x,Ncoef);
amin = min([a0; a]); amax = max([a0; a]);
bmin = min(b); bmax = max(b);
xg5 = fourier_series_gen(tp,T0,a0,a,b);

% Plot Coefficients
f = [f0:f0:(Ncoef*f0)]';
figure

% a's
%subplot(2,1,1),stem([0; f],[a0; a],'r','LineWidth',3,'MarkerSize',5), grid
%xlabel('f (Hz)');
%ylabel('a_m');
%title('Cosine Coefficients');
%axis([-f0 (Ncoef+1)*f0 amin amax]);
% b's
%subplot(2,1,2), stem(f,b,'LineWidth',3,'MarkerSize',5), grid
%xlabel('f (Hz)');
%ylabel('b_m');
%title('Sine Coefficients');
%axis([-f0 (Ncoef+1)*f0 bmin bmax]);

% 15
Ncoef = 15;
[a0,a,b]=fourier_series_coef(t,x,Ncoef);
xg15 = fourier_series_gen(tp,T0,a0,a,b);

% 50
Ncoef = 50;
[a0,a,b]=fourier_series_coef(t,x,Ncoef);
xg50 = fourier_series_gen(tp,T0,a0,a,b);
figure, plot(tp,xg5,tp,xg15,tp,xg50), grid on
xlabel('t (sec)');
ylabel('\it{x''(t)}')
title('Partial Sum Approximations (5, 15, 50 terms)');
legend('5','15','50')

%plotting Dn
%n=[-10:10];
subplot(2,1,1),stem([fftshift(Dn)],'r','LineWidth',3,'MarkerSize',5), grid
xlabel('n');
ylabel('magnitude');
title('magnitude of Dn');
axis([-f0 (Ncoef+1)*f0 amin amax]);

subplot(2,1,2), stem([angle(Dn(1:NNNd))],'LineWidth',3,'MarkerSize',5), grid
xlabel('n');
ylabel('argument');
title('argument of Dn');
axis([-f0 (Ncoef+1)*f0 bmin bmax]);

Function m-files:
fourier_series_coef.m

function [a0, a, b]=fourier_series_coef(t,x,numterms)
[nrt nct]=size(t);
[nrx ncx]=size(x);

% Check for input errors
if (nrt~=nrx)|(nct~=ncx)
error('t and x must be the same size')
end

if min(nrt,nct)~=1
error('t and x must be row or column vectors')
end
N = length(t);

td = diff(t);
if any(td<0)
error('t must be increasing')
end
tdd=diff(td);

dt = t(2)-t(1);
T0 = t(N)-t(1);
f0 = 1/T0;
w0 = 2*pi*f0;

% a0 = 1/T0 * integ(x(t)) over T0

a0 = (1/T0)*trapz(t,x);

% am = 2/T0 * integ(x(t)*cos(m*w0*t)) over T0
% bm = 2/T0 * integ(x(t)*sin(m*w0*t)) over T0

a = zeros(numterms,1);
b = zeros(numterms,1);
for m=1:numterms
a(m) = (2/T0)*trapz(t,(x.*cos(m*w0*t)));
b(m) = (2/T0)*trapz(t,(x.*sin(m*w0*t)));
end

return

fourier_series_gen.m

function [x] = fourier_series_gen(t,T0,a0,a,b)
% generates a waveform using the partial sum of the fourier series

[nra nca]=size(a);
[nrb ncb]=size(b);
Na = length(a);
Nb = length(b);

x = a0 + zeros(size(t));
w0=2*pi/T0;
for m=1:Na
x = x + a(m)*cos(m*w0*t);
end
for m=1:Nb
x = x + b(m)*sin(m*w0*t);
end

return

rect_fcn.m

function y = rect_fcn(t)
% Continuous-time rect function
% y = 0.0 for t<-0.5 and t>+0.5,
% y = 1.0 for -0.5< t <+0.5
% y = 0.5 for |t| = 0.5

y = zeros(size(t));
y((t>-0.5)&(t<+0.5)) = 1.0;
y((t==-0.5)|(t==+0.5)) = 0.5;
return