I am working on a project where I am trying to graph the output of input various functions represented by fourier series, applied to a transfer function of a bandreject filter.
I am getting my output waveform "flipped" from what I would expect it to be (see the graph of the sawtooth wave, but can't see what I've done wrong!
I noticed that by placing a negative sign in front of H on line 236 (the third line after "Sawtooth Wave Output" comment, "lambda_saw = H_saw.*C_n_saw;") it looks how I would expect, but I'm mystified as to why. Thank you for any help or comments, I greatly appreciate it.
%% Values (for analysis)
R = 1000/pi;
C = 100e-9;
sigma = 0.95;
%% Transfer Function (currently unused)
a = [(R*C)^2 4*R*C*(1-sigma) 1]; % denominator
b = [(R*C)^2 0 1]; % numerator
%% Part 3.1: Cosinusoids
%% Global
C_n_cos = zeros(1,201); % X corresponds with C_n
C_n_cos(1,100) = 1/2;
C_n_cos(1,102) = 1/2;
n = -100:100;
C_nM_cos = repmat(C_n_cos.',[1, 1e4]); % only needs to declared once - identical for all frequencies. 1e4 is size of t vector
%% 5Hz Input
T_0 = 1/5;
t= linspace(0, 3*T_0 ,1e4);
t_5Hz = t;
basis_5Hz = exp(1i*1e1*pi*n.'*t);
cos_5Hz = sum(C_nM_cos.*basis_5Hz);
%% 5 Hz Output
omega = (2*pi/T_0).*n; % omega is a vector
H_5Hz = freqs(b,a,omega);
lambda_5Hz = H_5Hz.*C_n_cos;
lambdaM_5Hz = repmat(lambda_5Hz',[1, length(t)]);
cos_5Hz_out = sum(lambdaM_5Hz.*basis_5Hz);
%% 50 Hz Input
T_0 = 1/50;
t= linspace(0, 3*T_0 ,1e4);
t_50Hz = t;
basis_50Hz = exp(1i*1e2*pi*n.'*t);
cos_50Hz = sum(C_nM_cos.*basis_50Hz);
%% 50 Hz Output
omega = (2*pi/T_0).*n; % omega is a vector
H_50Hz = freqs(b,a,omega);
lambda_50Hz = H_50Hz.*C_n_cos;
lambdaM_50Hz = repmat(lambda_50Hz',[1, length(t)]);
cos_50Hz_out = sum(lambdaM_50Hz.*basis_50Hz);
%% 500 Hz Input
T_0 = 1/500;
t= linspace(0, 3*T_0 ,1e4);
t_500Hz = t;
basis_500Hz = exp(1i*1e3*pi*n.'*t);
cos_500Hz = sum(C_nM_cos.*basis_500Hz);
%% 500 Hz Output
omega = (2*pi/T_0).*n; % omega is a vector
H_500Hz = freqs(b,a,omega);
lambda_500Hz = H_500Hz.*C_n_cos;
lambdaM_500Hz = repmat(lambda_500Hz',[1, length(t)]);
cos_500Hz_out = sum(lambdaM_500Hz.*basis_500Hz);
%% 5 kHz
T_0 = 1/(5*1e3);
t= linspace(0, 3*T_0 ,1e4);
t_5kHz = t;
basis_5kHz = exp(1i*1e4*pi*n.'*t);
cos_5kHz = sum(C_nM_cos.*basis_5kHz);
%% 5 kHz Output
omega = (2*pi/T_0).*n; % omega is a vector
H_5kHz = freqs(b,a,omega);
lambda_5kHz = H_5kHz.*C_n_cos;
lambdaM_5kHz = repmat(lambda_5kHz',[1, length(t)]);
cos_5kHz_out = sum(lambdaM_5kHz.*basis_5kHz);
%% 50 kHz
T_0 = 1/(5*1e4);
t= linspace(0, 3*T_0, 1e4);
t_50kHz = t;
basis_50kHz = exp(1i*100000*pi*n.'*t);
cos_50kHz = sum(C_nM_cos.*basis_50kHz);
%% 50 kHz Output
omega = (2*pi/T_0).*n; % omega is a vector
H_50kHz = freqs(b,a,omega);
lambda_50kHz = H_50kHz.*C_n_cos;
lambdaM_50kHz = repmat(lambda_50kHz',[1, length(t)]);
cos_50kHz_out = sum(lambdaM_50kHz.*basis_50kHz);
%% 500 kHz
T_0 = 1/(5*1e5);
t= linspace(0, 3*T_0 ,1e4);
t_500kHz = t;
basis_500kHz = exp(1i*1e6*pi*n.'*t);
cos_500kHz = sum(C_nM_cos.*basis_500kHz);
%% 500 kHz Output
omega = (2*pi/T_0).*n; % omega is a vector
H_500kHz = freqs(b,a,omega);
lambda_500kHz = H_500kHz.*C_n_cos;
lambdaM_500kHz = repmat(lambda_500kHz',[1, length(t)]);
cos_500kHz_out = sum(lambdaM_500kHz.*basis_500kHz);
%% 5 MHz Input
T_0 = 1/(5*1e6);
t= linspace(0, 3*T_0 ,1e4);
t_5MHz = t;
basis_5MHz = exp(1i*1e7*pi*n.'*t);
cos_5MHz = sum(C_nM_cos.*basis_5MHz);
%% 5 MHz Output
omega = (2*pi/T_0).*n; % omega is a vector
H_5MHz = freqs(b,a,omega);
lambda_5MHz = H_5MHz.*C_n_cos;
lambdaM_5MHz = repmat(lambda_5MHz',[1, length(t)]);
cos_5MHz_out = sum(lambdaM_5MHz.*basis_5MHz);
%% Cosinusoid Plots
figure('Name', '1 Volt Cosinusiods, 5 Hz to 5 MHz', 'NumberTitle','off');
subplot (2,4,1)
plot(t_5Hz, cos_5Hz, t_5Hz, cos_5Hz_out);
title('Cosinusoid, f = 5 Hz, T = 0.2 s');
xlabel('Time [s]'); ylabel('Voltage [V]');
subplot (2,4,2)
plot(t_50Hz, cos_50Hz, t_50Hz, cos_50Hz_out);
title('Cosinusoid, f = 50 Hz, T = 0.02 s');
xlabel('Time [s]'); ylabel('Voltage [V]');
subplot (2,4,3)
plot(t_500Hz, cos_500Hz, t_500Hz, cos_500Hz_out);
title('Cosinusoid, f = 500 Hz, T = 2 ms');
xlabel('Time [s]'); ylabel('Voltage [V]');
subplot (2,4,4)
plot(t_5kHz, cos_5kHz, t_5kHz, cos_5kHz_out);
title('Cosinusoid, f = 5 kHz, T = 0.2 ms');
xlabel('Time [s]'); ylabel('Voltage [V]');
subplot (2,4,5)
plot(t_50kHz, cos_50kHz, t_50kHz, cos_50kHz_out);
title('Cosinusoid, f = 50 kHz, T = 20 \mus');
xlabel('Time [s]'); ylabel('Voltage [V]');
subplot (2,4,6)
plot(t_500kHz, cos_500kHz, t_500kHz, cos_500kHz_out);
title('Cosinusoid, f = 500 kHz, T = 2 \mus');
xlabel('Time [s]'); ylabel('Voltage [V]');
subplot (2,4,7)
plot(t_5MHz, cos_5MHz, t_5MHz, cos_5MHz_out);
title('Cosinusoid, f = 5MHz, T = 200 ns');
xlabel('Time [s]'); ylabel('Voltage [V]');
%% Part 3.2: Square Wave
%% Square Wave Input
T_0 = 1/5000; %fundamental frequency, check.
t = linspace(0,T_0*3,1e4); % time vector [s]
n = -100:100; % indices included in the summation
C_n_square = 5./pi./n.*sin(pi*n/2); % this constructs the coefficients of
% e^(j*2*pi*n*t/T) for the sum; i.e., the c_k's
C_n_square(n==0) = 2.5; % note that the n = 0 is outside the
% summation in our representation
C_nM_square = repmat(C_n_square.',[1, length(t)]); % this replicates the transpose of the
% matrix X as many times as t is long
basis = exp(1i*10000*pi*n.'*t); % this is the set of exponentials in the sum.
% Note that this is a matrix as the value of
% the function depends *both* on n and t
SquareWaveIn = sum(C_nM_square.*basis); % 'sum' will sum down the columns of the matrix
% formed by element-wise multiplication of XM
% and basis
if max(imag(SquareWaveIn))<2e-15 % this checks to see if the imaginary part of x
SquareWaveIn = real(SquareWaveIn); % is below Matlab's noise floor (approx 10^-15).
end % If so, we get rid of the imaginary part.
%% Square Wave Output
omega = (2*pi/T_0).*n; % omega is a vector
H_square = freqs(b,a,omega);
lambda_square = H_square.*C_n_square;
lambdaM_square = repmat(lambda_square',[1, length(t)]);
SquareWaveOut = sum(lambdaM_square.*basis);
if max(imag(SquareWaveOut))<2e-15 % this checks to see if the imaginary part of x
SquareWaveOut = real(SquareWaveOut); % is below Matlab's noise floor (approx 10^-15).
end % If so, we get rid of the imaginary part.
%% Square Wave Plots
% we will use these variables to keep the limits on our plots uniform
tmin = min(t); tmax = max(t); xmin = min(SquareWaveOut - 1); xmax = max(SquareWaveOut + 1);
figure('Name', 'Square Wave, N = 100 Fourier Approximation', 'NumberTitle','off');
plot(t,SquareWaveIn,t,SquareWaveOut);
axis([tmin tmax xmin xmax]); % sets the plot's axis limits
title('Square Wave, f = 5 kHz, T = 200 \mus');
xlabel('Time [s]'); ylabel('Voltage [V]');
%% Part 3.3 : Sawtooth Wave
T_0 = 1/2500;
t = linspace(0,T_0*3,1e4); % time vector [s]
n = -100:100; % indices included in the summation, n is a 101 element vector
%C_n_saw = (5*1i)./(pi*n); % this constructs the coefficients of e^(*) terms
C_n_saw = (1i*5*((-1).^n))./(pi.*n);
C_n_saw(n==0) = 0; % C_0
% X is a horizontal vector with the coefficients of the exponential term
% corresponding each n-value from -100 to 100.
C_nM_saw = repmat(C_n_saw.',[1, length(t)]); % This is a matrix that replicates the transpose of the
% vector X as many times as t is long
% now, XM is a matrix with identical values in each row.
% Each column has a coefficient of an exponential term corresponding to an
% n-value, repeated throughout
basis = exp(1i*pi*5000*n.'*t); % basis is a matrix containing the exponential part,
% each column corresponding with each n-value
% and each row corresponding time
SawtoothWaveIn = sum(C_nM_saw.*basis); % x100 is matrix resulting from summing down the rows of element-wise
% multiplacation of XM and basis
if max(imag(SawtoothWaveIn))<3e-15 % this checks to see if the imaginary part of x
SawtoothWaveIn = real(SawtoothWaveIn); % is below Matlab's noise floor (approx 10^-15).
end % If so, we get rid of the imaginary part.
%% Sawtooth Wave Output
omega = (2*pi/T_0).*n; % omega is a vector
H_saw = freqs(b,a,omega);
lambda_saw = H_saw.*C_n_saw;
lambdaM_saw = repmat(lambda_saw',[1, length(t)]);
SawtoothWaveOut = sum(lambdaM_saw.*basis);
if max(imag(SawtoothWaveOut))<3e-15 % this checks to see if the imaginary part of x
SawtoothWaveOut = real(SawtoothWaveOut); % is below Matlab's noise floor (approx 10^-15).
end % If so, we get rid of the imaginary part.
%% Sawtooth Wave Plots
% we will use these variables to keep the limits on our plots uniform
tmin = min(t); tmax = max(t); xmin = min(SawtoothWaveOut - 1); xmax = max(SawtoothWaveOut + 1);
figure('Name', 'Sawtooth Wave, n = 100 Fourier Approximation', 'NumberTitle','off');
plot(t, SawtoothWaveIn, t, SawtoothWaveOut);
axis([tmin tmax xmin xmax]); % sets the plot's axis limits
title('Sawtooth Wave, f = 2.5 kHz, T = 400 \mus');
xlabel('Time [s]'); ylabel('Voltage [V]');
C_n_sawneeds a different expression (even though its FFT shows a ramp)? Maybe it derives from exp(-ix), rather than exp(ix)? \$\endgroup\$