Assignment:
Part 1:
Step 1: Download the monthly sunspot data for the 33-year period from January 1875 through December 1907. This will be a total of 396 data points.
Go to https://sidc.oma.be/sunspot-data/.
Under ASCII files, select monthly and monthly smoothed sunspot number . Save the data to the current directory of your MATLAB work space (denoted at the top when you open MATLAB) as a text file.
Open the file by double-clicking and adjust the Notepad window until four columns display. Then delete everything from the top down to where the first 1875 appears. Delete everything from after the last 1907 down to the bottom of the file. Save the resulting file.
Import the data into MATLAB, type the following at the MATLAB prompt:
>> [a b c d e]=textread('monthssn.txt','%f %f %f %c %u');
If your file is not saved in your current directory, you will have to type the complete path. The data we are interested in is in array b. We can rename this spd by typing the following at the prompt:
>>spd=b;
If you plot spd, you will see that the data does indeed contain a cyclic component.
Step 2: Remove the average value of spd from the data so that the constant component in the Fourier transform is removed.
Question 1: What is the frequency of the cyclic component in the data.?
Question 2: How can we be sure that the average value was properly removed from the data?
Part 2:
Step 1: Generate random noise of mean (average) zero and 1024 points:
x=rand(1024,1)-0.5;
The command rand generates random numbers uniformly distributed between 0 and 1. Thus, we subtract 0.5 to obtain a mean of zero.
Step 2: Compute and plot the amplitude spectrum for this discrete signal.
Question 1: What frequencies should we expect to see in the spectrum?
Part 3:
(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.
(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(w) and an approximation of Y to Y(w). Plot |X(w)| and |Y(w)| versus w , 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. We only need the frequencies from 0 to the Nyquist frequency, so only consider the plot in the range of 0 ≤w≤30 rad/sec.
(d). Recall that Y(w) = X(w)H(w). Compute the approximation to H(w) by dividing the elements of Y by the elements of X. Plot the frequency response |H(w)|.
Questions: What sort of filter does this represent? What is its bandwidth?
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 nx = 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;