SHARE
TWEET

Upsampling Interpolation - various filters

a guest Apr 29th, 2018 141 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. pkg load signal
  2.  
  3. % big values take too long for unoptimized calculation, so
  4. % for illustration, we will use 44.1 instead of 44100 (divide all by 1000)
  5.  
  6. Fs = 44.1; % the sampling rate
  7. upsamplingFactor = 16; % 16X upsamp+interp. This matches Blu2 44.1kHz to 705.6kHz
  8. hz = 19; % the sine wave will be at this frequency
  9.          % this must be < Fs/2 to avoid aliasing
  10. hz1 = hz; % machinery actually can start at hz1 and end at hz2 for more illustrations...
  11. hz2 = hz;
  12.  
  13.          
  14. F = upsamplingFactor*Fs; % the upsampling frequency
  15. Fdisp = min(128,upsamplingFactor*4)*Fs; % resolution of the continuous plot
  16. T = 1/Fs; % the sampling time
  17.  
  18. tmax = 5;  % generate and calculate from 0 to this time
  19. tdisp = 0:1/Fdisp:tmax; % the time-samples for "continuous" plotting
  20. t = 0:1/F:tmax; % the time-samples for upsampled plotting
  21. ts = 0:T:tmax;  % the time at the sample points
  22.  
  23. % basic sine
  24. %x = @(t) sin(2*pi*hz*t);
  25.  
  26. % sine sweep -- frequency varies from hz1 at t=0
  27. %     to hz2 at t=tmax
  28. x = @(t) sin(2*pi*(t.*t/tmax*hz2 + (t-t.*t/tmax)*hz1));
  29.  
  30. xn = x(ts);  % the discrete sequence
  31. xf = x(t);
  32.  
  33. tmiddle = 0.2; % plot this much time from the center
  34.  
  35. width = 2;
  36. height = 6;
  37. num = 1;
  38. fig = 1;
  39. lpfFreq = 22;
  40.  
  41. % the wave with the samples
  42. subplot(height,width,num); num=num+1;
  43. hold off;
  44. plot(tdisp, x(tdisp), 'k'); hold on;  % plot the sine wave
  45. stem(ts, xn);  hold on; % plot sampled
  46. xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
  47.                                        %   (bec the edges are not bw-limited)
  48. ylim([-1,1]);
  49. title(strcat(num2str(fig),'. Original Wave with Sample Points -',{' '},
  50.       num2str(hz),' Hz'));
  51.  
  52. % just the samples
  53. subplot(height,width,num); num=num+1;
  54. hold off;
  55. stem(ts, xn);
  56. xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
  57.                                        %   (bec the edges are not bw-limited)
  58. ylim([-1,1]);
  59. fig = fig + 1;
  60. title(strcat(num2str(fig),'. Sampled Points Only - Sampled at',{' '},
  61.       num2str(Fs),' Hz'));
  62.  
  63. num = num+1;
  64. % perform interpolation using sinc:
  65. interpolated = 0;
  66. for n=1:length(ts)
  67.     % n-1 due to the 1-based indexing of matlab
  68.     interpolated = interpolated + xn(n) * sinc((t-(n-1)*T)/T);
  69. end
  70. %%sse = sum((interpolated-xf)(:).^2);
  71.  
  72. subplot(height,width,num); num=num+1;
  73. hold off;
  74. stem(ts, xn);  hold on; % plot sampled
  75. stem(t, interpolated, 'k:.'); % plot interpolated in black
  76. xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
  77.                                        %   (bec the edges are not bw-limited)
  78. ylim([-1,1]);
  79. fig = fig+1;
  80. line1 = strcat(num2str(fig),'. Whittaker-Shannon - Ideal (Infinite) Sinc to get',
  81.   {' '},num2str(upsamplingFactor),'X upsamp+interp');
  82. title([line1,'(Ideal brickwall filter, i.e. perfect rectangle in freq domain)']);
  83.  
  84.  
  85. % using zero order hold (that is, repeating the samples)
  86. %  (truncated end to match t length)
  87. %interpolated =  repmat(xn(:)', F/Fs)(1:(length(xn)-1)*F/Fs+1);
  88. interpolated =  repmat(xn(:)', F/Fs)(1:length(t));
  89. subplot(height,width,num); num=num+1;
  90. hold off;
  91. stem(ts, xn);  hold on; % plot sampled
  92. stem(t, interpolated, 'k:.');
  93. xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
  94.                                        %   (bec the edges are not bw-limited)
  95. ylim([-1,1]);
  96. fig = fig+1;
  97. line1 = strcat(num2str(fig),'. Repeating Each Sample to get',{' '},
  98.       num2str(upsamplingFactor),'X upsamp+interp');
  99. title([line1,'(Zero-order hold)']);
  100.  
  101.  
  102. % Connect-the-dots: Draw a line between 2 adjacent samples
  103. interpolated =  interp1(ts,xn,t);
  104. adjEnd = length(xf)-upsamplingFactor; % interp1 has NA at the end bec there's
  105.                                       %   no next sample to draw a line to
  106. subplot(height,width,num); num=num+1;
  107. hold off;
  108. stem(ts, xn);  hold on; % plot sampled
  109. stem(t, interpolated, 'k:.');
  110. xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
  111.                                        %   (bec the edges are not bw-limited)
  112. ylim([-1,1]);
  113. fig = fig+1;
  114. title(strcat(num2str(fig),'. Connect the Dots to get',{' '},
  115.       num2str(upsamplingFactor),'X upsamp+interp'));
  116.  
  117.  
  118. % simple FIR LPF
  119. lpfOrder = 100;
  120. f =  [lpfFreq ]/(F/2);
  121. hc = fir1(lpfOrder, f,'low');
  122. zeroStuffed = upsample(xn,upsamplingFactor);
  123. interpolated = upsamplingFactor*filter(hc,1,zeroStuffed);
  124. delay = floor(lpfOrder / 2); % how do we shift the signal.... ts no longer lines up
  125. interpolated = interpolated(1:length(t));
  126.  
  127. subplot(height,width,num); num=num+1;
  128. hold off;
  129. plot((-0.5:1/4096:0.5-1/4096)*F,20*log10(abs(fftshift(fft(hc,4096)))));
  130. hold on;
  131. plot((-0.5:1/4096:0.5-1/4096)*F,
  132.       20*log10(abs(fftshift(fft(zeroStuffed,4096))))-80);
  133. xlim([0,Fs*2]);
  134. ylim([-100,1]);
  135. fig = fig + 1;
  136. line1 = strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
  137.     ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter Frequency Response in blue;');
  138. title([line1,
  139.     'upsampled (but not interpolated) signal in orange']);
  140. grid on;
  141.  
  142. subplot(height,width,num); num=num+1;
  143. hold off;
  144. stem(ts, xn);  hold on; % plot sampled
  145. hold on;
  146. stem(t-delay/F, interpolated, 'k:.'); % plot interpolated in black
  147. xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
  148.                                        %   (bec the edges are not bw-limited)
  149. ylim([-1,1]);
  150. fig = fig + 1;
  151. title(strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
  152.     ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter to get',
  153.     {' '},num2str(upsamplingFactor),'X upsamp+interp'));
  154.    
  155.  
  156.  
  157. % simple FIR LPF 2
  158. lpfOrder = 164;
  159. f =  [lpfFreq ]/(F/2);
  160. hc = fir1(lpfOrder, f,'low');
  161. zeroStuffed = upsample(xn,upsamplingFactor);
  162. interpolated = upsamplingFactor*filter(hc,1,zeroStuffed);
  163. delay = floor(lpfOrder / 2); % how do we shift the signal.... ts no longer lines up
  164. interpolated = interpolated(1:length(t));
  165.  
  166. subplot(height,width,num); num=num+1;
  167. hold off;
  168. plot((-0.5:1/4096:0.5-1/4096)*F,20*log10(abs(fftshift(fft(hc,4096)))));
  169. hold on;
  170. plot((-0.5:1/4096:0.5-1/4096)*F,
  171.       20*log10(abs(fftshift(fft(zeroStuffed,4096))))-80);
  172. xlim([0,Fs*2]);
  173. ylim([-100,1]);
  174. fig = fig + 1;
  175. line1 = strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
  176.     ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter Frequency Response in blue;');
  177. title([line1,
  178.     'upsampled (but not interpolated) signal in orange']);
  179. grid on;
  180.  
  181. subplot(height,width,num); num=num+1;
  182. hold off;
  183. stem(ts, xn);  hold on; % plot sampled
  184. hold on;
  185. stem(t-delay/F, interpolated, 'k:.'); % plot interpolated in black
  186. xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
  187.                                        %   (bec the edges are not bw-limited)
  188. ylim([-1,1]);
  189. fig = fig + 1;
  190. title(strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
  191.     ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter to get',
  192.     {' '},num2str(upsamplingFactor),'X upsamp+interp'));
  193.    
  194.  
  195. % simple FIR LPF 3
  196. lpfOrder = 1016;
  197. f =  [lpfFreq ]/(F/2);
  198. hc = fir1(lpfOrder, f,'low');
  199. zeroStuffed = upsample(xn,upsamplingFactor);
  200. interpolated = upsamplingFactor*filter(hc,1,zeroStuffed);
  201. delay = floor(lpfOrder / 2); % how do we shift the signal.... ts no longer lines up
  202. interpolated = interpolated(1:length(t));
  203.  
  204. subplot(height,width,num); num=num+1;
  205. hold off;
  206. plot((-0.5:1/4096:0.5-1/4096)*F,20*log10(abs(fftshift(fft(hc,4096)))));
  207. hold on;
  208. plot((-0.5:1/4096:0.5-1/4096)*F,
  209.       20*log10(abs(fftshift(fft(zeroStuffed,4096))))-80);
  210. xlim([0,Fs*2]);
  211. ylim([-100,1]);
  212. fig = fig + 1;
  213. line1 = strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
  214.     ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter Frequency Response in blue;');
  215. title([line1,
  216.     'upsampled (but not interpolated) signal in orange']);
  217. grid on;
  218.  
  219. subplot(height,width,num); num=num+1;
  220. hold off;
  221. stem(ts, xn);  hold on; % plot sampled
  222. hold on;
  223. stem(t-delay/F, interpolated, 'k:.'); % plot interpolated in black
  224. xlim([tmax/2-tmiddle/2,tmax/2+tmiddle/2]); % pick graph from the middle
  225.                                        %   (bec the edges are not bw-limited)
  226. ylim([-1,1]);
  227. fig = fig + 1;
  228. title(strcat(num2str(fig),'. FIR ',{' '},num2str(lpfOrder),
  229.     ' order,',{' '},num2str(lpfFreq),'Hz Low Pass Filter to get',
  230.     {' '},num2str(upsamplingFactor),'X upsamp+interp'));
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top