Upsampling Interpolation - various filters

pkg load signal

% big values take too long for unoptimized calculation, so
% for illustration, we will use 44.1 instead of 44100 (divide all by 1000)

Fs = 44.1; % the sampling rate
upsamplingFactor = 16; % 16X upsamp+interp. This matches Blu2 44.1kHz to 705.6kHz
hz = 19; % the sine wave will be at this frequency
         % this must be < Fs/2 to avoid aliasing
hz1 = hz; % machinery actually can start at hz1 and end at hz2 for more illustrations...
hz2 = hz;


F = upsamplingFactor*Fs; % the upsampling frequency
Fdisp = min(128,upsamplingFactor*4)*Fs; % resolution of the continuous plot
T = 1/Fs; % the sampling time

tmax = 5;  % generate and calculate from 0 to this time
tdisp = 0:1/Fdisp:tmax; % the time-samples for "continuous" plotting
t = 0:1/F:tmax; % the time-samples for upsampled plotting
ts = 0:T:tmax;  % the time at the sample points

% basic sine
%x = @(t) sin(2*pi*hz*t);

% sine sweep -- frequency varies from hz1 at t=0
%     to hz2 at t=tmax
x = @(t) sin(2*pi*(t.*t/tmax*hz2 + (t-t.*t/tmax)*hz1));

xn = x(ts);  % the discrete sequence
xf = x(t);

tmiddle = 0.2; % plot this much time from the center

width = 2;
height = 6;
num = 1;
fig = 1;
lpfFreq = 22;

% the wave with the samples
subplot(height,width,num); num=num+1;
hold off;
plot(tdisp, x(tdisp), 'k'); hold on;  % plot the sine wave
stem(ts, xn);  hold on; % plot sampled
xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
                                       %   (bec the edges are not bw-limited)
ylim([-1,1]);
title(strcat(num2str(fig),'. Original Wave with Sample Points -',{' '},
      num2str(hz),' Hz'));

% just the samples
subplot(height,width,num); num=num+1;
hold off;
stem(ts, xn);
xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
                                       %   (bec the edges are not bw-limited)
ylim([-1,1]);
fig = fig + 1;
title(strcat(num2str(fig),'. Sampled Points Only - Sampled at',{' '},
      num2str(Fs),' Hz'));

num = num+1;
% perform interpolation using sinc:
interpolated = 0;
for n=1:length(ts)
    % n-1 due to the 1-based indexing of matlab
    interpolated = interpolated + xn(n) * sinc((t-(n-1)*T)/T);
end
%%sse = sum((interpolated-xf)(:).^2);

subplot(height,width,num); num=num+1;
hold off;
stem(ts, xn);  hold on; % plot sampled
stem(t, interpolated, 'k:.'); % plot interpolated in black
xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
                                       %   (bec the edges are not bw-limited)
ylim([-1,1]);
fig = fig+1;
line1 = strcat(num2str(fig),'. Whittaker-Shannon - Ideal (Infinite) Sinc to get',
  {' '},num2str(upsamplingFactor),'X upsamp+interp');
title([line1,'(Ideal brickwall filter, i.e. perfect rectangle in freq domain)']);


% using zero order hold (that is, repeating the samples)
%  (truncated end to match t length)
%interpolated =  repmat(xn(:)', F/Fs)(1:(length(xn)-1)*F/Fs+1);
interpolated =  repmat(xn(:)', F/Fs)(1:length(t));
subplot(height,width,num); num=num+1;
hold off;
stem(ts, xn);  hold on; % plot sampled
stem(t, interpolated, 'k:.');
xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
                                       %   (bec the edges are not bw-limited)
ylim([-1,1]);
fig = fig+1;
line1 = strcat(num2str(fig),'. Repeating Each Sample to get',{' '},
      num2str(upsamplingFactor),'X upsamp+interp');
title([line1,'(Zero-order hold)']);


% Connect-the-dots: Draw a line between 2 adjacent samples
interpolated =  interp1(ts,xn,t);
adjEnd = length(xf)-upsamplingFactor; % interp1 has NA at the end bec there's
                                      %   no next sample to draw a line to
subplot(height,width,num); num=num+1;
hold off;
stem(ts, xn);  hold on; % plot sampled
stem(t, interpolated, 'k:.');
xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
                                       %   (bec the edges are not bw-limited)
ylim([-1,1]);
fig = fig+1;
title(strcat(num2str(fig),'. Connect the Dots to get',{' '},
      num2str(upsamplingFactor),'X upsamp+interp'));


% simple FIR LPF
lpfOrder = 100;
f =  [lpfFreq ]/(F/2);
hc = fir1(lpfOrder, f,'low');
zeroStuffed = upsample(xn,upsamplingFactor);
interpolated = upsamplingFactor*filter(hc,1,zeroStuffed);
delay = floor(lpfOrder / 2); % how do we shift the signal.... ts no longer lines up
interpolated = interpolated(1:length(t));

subplot(height,width,num); num=num+1;
hold off;
plot((-0.5:1/4096:0.5-1/4096)*F,20*log10(abs(fftshift(fft(hc,4096)))));
hold on;
plot((-0.5:1/4096:0.5-1/4096)*F,
      20*log10(abs(fftshift(fft(zeroStuffed,4096))))-80);
xlim([0,Fs*2]);
ylim([-100,1]);
fig = fig + 1;
line1 = strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
    ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter Frequency Response in blue;');
title([line1,
    'upsampled (but not interpolated) signal in orange']);
grid on;

subplot(height,width,num); num=num+1;
hold off;
stem(ts, xn);  hold on; % plot sampled
hold on;
stem(t-delay/F, interpolated, 'k:.'); % plot interpolated in black
xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
                                       %   (bec the edges are not bw-limited)
ylim([-1,1]);
fig = fig + 1;
title(strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
    ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter to get',
    {' '},num2str(upsamplingFactor),'X upsamp+interp'));



% simple FIR LPF 2
lpfOrder = 164;
f =  [lpfFreq ]/(F/2);
hc = fir1(lpfOrder, f,'low');
zeroStuffed = upsample(xn,upsamplingFactor);
interpolated = upsamplingFactor*filter(hc,1,zeroStuffed);
delay = floor(lpfOrder / 2); % how do we shift the signal.... ts no longer lines up
interpolated = interpolated(1:length(t));

subplot(height,width,num); num=num+1;
hold off;
plot((-0.5:1/4096:0.5-1/4096)*F,20*log10(abs(fftshift(fft(hc,4096)))));
hold on;
plot((-0.5:1/4096:0.5-1/4096)*F,
      20*log10(abs(fftshift(fft(zeroStuffed,4096))))-80);
xlim([0,Fs*2]);
ylim([-100,1]);
fig = fig + 1;
line1 = strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
    ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter Frequency Response in blue;');
title([line1,
    'upsampled (but not interpolated) signal in orange']);
grid on;

subplot(height,width,num); num=num+1;
hold off;
stem(ts, xn);  hold on; % plot sampled
hold on;
stem(t-delay/F, interpolated, 'k:.'); % plot interpolated in black
xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
                                       %   (bec the edges are not bw-limited)
ylim([-1,1]);
fig = fig + 1;
title(strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
    ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter to get',
    {' '},num2str(upsamplingFactor),'X upsamp+interp'));


% simple FIR LPF 3
lpfOrder = 1016;
f =  [lpfFreq ]/(F/2);
hc = fir1(lpfOrder, f,'low');
zeroStuffed = upsample(xn,upsamplingFactor);
interpolated = upsamplingFactor*filter(hc,1,zeroStuffed);
delay = floor(lpfOrder / 2); % how do we shift the signal.... ts no longer lines up
interpolated = interpolated(1:length(t));

subplot(height,width,num); num=num+1;
hold off;
plot((-0.5:1/4096:0.5-1/4096)*F,20*log10(abs(fftshift(fft(hc,4096)))));
hold on;
plot((-0.5:1/4096:0.5-1/4096)*F,
      20*log10(abs(fftshift(fft(zeroStuffed,4096))))-80);
xlim([0,Fs*2]);
ylim([-100,1]);
fig = fig + 1;
line1 = strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
    ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter Frequency Response in blue;');
title([line1,
    'upsampled (but not interpolated) signal in orange']);
grid on;

subplot(height,width,num); num=num+1;
hold off;
stem(ts, xn);  hold on; % plot sampled
hold on;
stem(t-delay/F, interpolated, 'k:.'); % plot interpolated in black
xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
                                       %   (bec the edges are not bw-limited)
ylim([-1,1]);
fig = fig + 1;
title(strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
    ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter to get',
    {' '},num2str(upsamplingFactor),'X upsamp+interp'));