%% Set some parameters
L = 18;
NFFT = 2^18;
freq = [0:NFFT-1]/NFFT - 0.5;
n = [-NFFT/2:NFFT/2-1];

%% Make a signal
f0 = 3/16; f0 = round(f0*NFFT)/NFFT;
X = zeros(NFFT,1);
X(1:f0*NFFT) = [f0*NFFT:-1:1]/(f0*NFFT);
X(end-f0*NFFT+1:end) = ones(f0*NFFT,1);
x = fftshift(ifft(X));
figure(1); clf;
subplot(211);
stem(n,[real(x),imag(x)]);
xlim([-40 40]);
subplot(212);
plot(freq,fftshift(X));
axis([-0.5 0.5 -0.1 1.1]);
grid on;
pause;

%% Upsample the signal
z = x;
nz = [-length(z)/2:length(z)/2-1];
Z = abs(fftshift(fft(z,NFFT)));

U = 3;
K = NFFT*U;
f = [0:K-1]/K - 0.5;
a = zeros(length(z)*U,1);
a(1:U:end) = z;
na = [-length(a)/2:length(a)/2-1];
A = abs(fftshift(fft(a,K)));

fp = 0.9/(2*U);
fs = 1.1/(2*U);
f1 = 0.5*(fs+fp);
f2 = 0.5*(fs-fp);
m = [-100:100];
h = (1/U)*sinc(2*f1*m).*sinc(2*f2*m);
H = abs(fftshift(fft(h,K)));

b = conv(a,h,'same');
B = abs(fftshift(fft(b,K)));

figure(1); clf;
subplot(221);
stem(nz,[real(z),imag(z)]);
xlim([-40 40]);
subplot(222);
stem(na,[real(a),imag(a)]);
xlim([-40 40]);
subplot(223);
plot(freq,Z);
grid on;
subplot(224);
plot(f,A);
grid on;
pause;

%% Now low pass filter
figure(1); clf;
subplot(221);
stem(na,[real(a),imag(a)]);
xlim([-40 40]);
subplot(222);
stem(na,[real(b),imag(b)]);
xlim([-40 40]);
subplot(223);
plot(f,H,'r'); hold on;
plot(f,A); hold off;
grid on;
subplot(224);
plot(f,B);
grid on;
pause;

%% Downsample the signal
x = b;
n = [-length(x)/2:length(x)/2-1];
NFFT = 2^nextpow2(length(b));
freq = [0:NFFT-1]/NFFT - 0.5;

N = 8;
p = zeros(size(x));
p(1:N:end) = 1;
y = p.*x;         % Create zeros
z = y(1:N:end);   % Remove zeros
nz = n(1:N:end)/N;
Y = abs(fftshift(fft(y,NFFT)));
Z = abs(fftshift(fft(z,NFFT)));
figure(1); clf;
subplot(221);
plot(n,[real(x),imag(x)]); hold on;
stem(n,[real(y),imag(y)]); hold off;
xlim([-40 40]);
subplot(222);
stem(nz,[real(z),imag(z)]);
xlim([-40 40]);
subplot(223);
plot(freq,Y);
grid on;
subplot(224);
plot(freq,Z);
grid on;

