﻿ # Upsampling Interpolation - various filters

a guest
Apr 29th, 2018
189
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
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