Estoy trabajando en un proyecto en el que estoy tratando de graficar la salida de varias funciones de entrada representadas por series de Fourier, aplicadas a una función de transferencia de un filtro de banda.
Estoy haciendo que mi forma de onda de salida se "desvíe" de lo que esperaba que fuera (vea la gráfica de la onda de diente de sierra, ¡pero no puedo ver lo que hice mal!
Me di cuenta de que al colocar un signo negativo delante de H en la línea 236 (la tercera línea después del comentario de "Sawtooth Wave Output", "lambda_saw = H_saw. * C_n_saw;") se ve como lo esperaría, pero Estoy desconcertado de por qué. Gracias por cualquier ayuda o comentario, lo aprecio mucho.
%% 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]');