Mark's Homework 15

Problem

Make a MATLAB script to do four times oversampling and a filter so as to eliminate as much as possible the effect of the D/A converter that follows the interpolation filter.

MATLAB Script

The following is currently NON-WORKING. A working script will follow.

clear all;
close all;
N = 50; %Number of sample points for our original signal
tmax=2; %Time, in seconds, to extend our original signal. We then sample this signal N times
oversampling = 4; %Number of oversampling points
T=tmax/N; %Calculate T, the sample period
t=0:T:(N)*T;
x=sin(2*pi*t)+sin(2*pi*2*t)+sin(2*pi*5*t); %This is our signal to over sample
N1=(oversampling-1)*N;
X=fft(x);
X1 = (N1+N)/N*[X(1:N/2), zeros(1, N1), X(N/2+1:N)]; %This is our interpolating filter

x2= ifft(X1);

figure(7)
t1 = 0:N/(N+N1)*T:(N-1/(N1+N))*T;
plot(t1, x2, 'bo');
hold on
t=0:T:(N)*T;
plot(t, x, 'ro', t, x, 'r-');
hold off

f = -1/(2*T):1/(N*T):1/(2*T);

figure(8)
plot(f, abs(fftshift(x)))

%This is to cancel out the effects of the D/A Converter
sT = T/100; %sT is hwo many times we are going to sample our pulse p(t)
stmax = oversampling*T; %This is how wide our pulse has to be when we oversample
sN = stmax/sT; %This is the total number of sampling points for sampling p(t)
t = 0:sT:stmax-sT;
p = u(t+1e-9) - u(t-T/(2*oversampling)) + u(t - (oversampling*T-T/(2*oversampling))); %This is p(t) sampled
P = fft(p);
iP = 1./P; %Take 1/P(f) so we cancel the effects of the D/A converter (which is P(f))
numpoints = 2*T/sT;
%We can only go out to 1/(2*T) on the 1/P(f). The rest we need to zero out.
iP = [iP(1:numpoints) zeros(1, sN-numpoints*2-1) iP(sN-numpoints:sN)]
iP = iP(1:sN/(N1+N):sN); %Down sample 1/P(f) so we can element multiply with our interpolating filter

%This plots P(f)
f = -1/(2*sT):1/(sN*sT):1/(2*sT)-1/(sN*sT);
Ps = fftshift(P); %Shift P(f)
%Figure 4 plots P(f) shifted
figure(4)
plot(f, abs(Ps));
%Figure 5 plots P(f) shifted and zoomed in the f axis
figure(5)
plot(f(sN/2-numpoints+1:sN/2+numpoints), abs(Ps(sN/2-numpoints+1:sN/2+numpoints)));

X1 = X1 .* iP; %Element multiply our interpolating filter with our D/A converter effect canceller filter
x1 = ifft(X1);

%Figure 1 plots the original signal and our interpolated signal (with D/A
%converter effects cancelled.)
figure(1)
t1 = 0:N/(N+N1)*T:(N-1/(N1+N))*T;
plot(t1, x2, 'bo');
hold on
t=0:T:(N)*T;
plot(t, x, 'ro', t, x, 'r-');
hold off

%Figure 2 plots just our interpolated signal, but this time with lines
%connecting each point.
figure(2)
plot(t1, x1, 'b-');

%Figure 3 plots the frequency response of our original sampled signal.
figure(3)
Xs((N+1)/2:N) = X(1:(N+1)/2); %The fftshift didn't work for me
Xs(1:(N+1)/2) = X((N+1)/2:N); %So this manually shifts the FFT
f = -1/(2*T):1/(N*T):1/(2*T)-1/(N*T);
plot(f,abs(Xs))
%figure(4)
%plot(f, abs(fftshift(X1)))