Fourier Analysis of system in MatLab

Use the contfft program to identify the frequency response of a system from its inputs and outputs. The program contfft is shown below.Save it as contfft.m in your current
MATLAB directory as denoted at the top of the MATLAB window.

(a). Generate a random input signal x(t) in MATLAB by using the command rand. In particular, let x = rand(1000,1)-0.5. This will create a 1000 point vector of random numbers that represents the samples of x(t). Assume that the sampling period was T = 0.1, and create a corresponding time vector.

(b). Use MATLAB to determine a corresponding output y(t) for a particular system. In particular, type the command y=conv(exp(-t),x); to generate y(t). The array y will be 1999 elements long by virtue of how the convolution works. Use only the first 1000 elements. Plot x versus t and y versus t, and compare the two signals in terms of frequency content. (In the time domain, we characterize high-frequency content qualitatively by sharp transitions or quick motion of the signal.)

(c). Suppose that you have a spectrum analyzer available to identify the system using only the signals x(t) and y(t). Use contfft to compute the approximation X to X(ω) and an approximation of Y to Y(ω). Plot |X(ω)| and |Y(ω)| versus ω, and compare the difference in frequency content of these two signals. The frequency vector provided by the contfft command ranges from 0 to 62.769 rad/sec. However, since the sampling period is 0.1 sec, the Nyquist frequency is pi/T, 31.4159 rad/sec. Reconstruction cannot occur well for frequencies higher than this frequency. The frequencies above the Nyquist frequency can, on the other hand, be viewed as negative frequencies. A command fftshift can be used which reorders the output of the fft command such that the frequencies run from -pi/T, the negative Nyquist frequency, to pi/T, the Nyquist frequency. However, you then have to change the frequency vector for plotting on your own. For this exercise, we only need the frequencies from 0 to the Nyquist frequency, so only consider the plot in the range of 0 ≤ ω ≤ 30 rad/sec. (d). Recall that Y(ω) = X(ω)H(ω). Compute the approximation to H(ω) by dividing the elements of Y by the elements of X. Plot the frequency response |H(ω)|. What sort of filter does this represent? What is its bandwidth?

See Contifft program below.

Contfft Program

function [X,w] = contfft(x,T);
% CONTFFT [X,w] = contfft(x,T)
%
% Computes the Fourier transform of a continuous time signal
% using the FFT. The input is the sampled continuous
% time signal x and the sampling time T. The output is
% the Fourier transform X(w) and the frequency vector w.
%
[n,m] = size(x);
if n<m,
x = x';
end
Xn = fft(x);
N = length(x);
n = 0:N-1;
n(1) = eps;
X = (1-exp(-j*2*pi*n/N))./(j*2*pi*n/N/T).*Xn.';
w = 2*pi*n/N/T;