• API
• FAQ
• Tools
• Archive
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!
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.

Top