-
Notifications
You must be signed in to change notification settings - Fork 0
/
resamp.py
84 lines (77 loc) · 3.11 KB
/
resamp.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
#!/usr/bin/python3
# Copyright © 2019 Bart Massey
# [This program is licensed under the "MIT License"]
# Please see the file LICENSE in the source
# distribution of this software for license terms.
#
# Python reimplementation of http://www.nicholson.com/rhn/dsp.html#3
# BSD Licensed per author.
# Please see comment at end of file for original source and
# licensing information.
from math import pi, sin, cos
def resamp(x, indat, fmax, fsr, wnwdth):
"""
QDSS Windowed-Sinc ReSampling
x = new sample point location (relative to old indexes)
(e.g. every other integer for 0.5x decimation)
indat = original data array
fmax = low pass filter cutoff frequency
fsr = sample rate
wnwdth = width of windowed Sinc used as the low pass filter
Returns a filtered new sample point.
"""
alim = len(indat)
r_g = 2 * fmax / fsr # Calc gain correction factor
r_y = 0
for i in range(-wnwdth//2, wnwdth//2-1):
j = int(x + i) # Calc input sample index
# calculate von Hann Window. Scale and calculate Sinc
r_w = 0.5 - 0.5 * cos(2 * pi * (0.5 + (j - x) / wnwdth))
r_a = 2 * pi * (j - x) * fmax / fsr
r_snc = 1
if r_a != 0:
r_snc = sin(r_a) / r_a
if j >= 0 and j < alim:
r_y = r_y + r_g * r_w * r_snc * indat[j]
return r_y
# rem - QDSS Windowed-Sinc ReSampling subroutine in Basic
# rem
# rem - This function can also be used for interpolation of FFT results
# rem
# rem function parameters
# rem : x = new sample point location (relative to old indexes)
# rem (e.g. every other integer for 0.5x decimation)
# rem : indat = original data array
# rem : alim = size of data array
# rem : fmax = low pass filter cutoff frequency
# rem : fsr = sample rate
# rem : wnwdth = width of windowed Sinc used as the low pass filter
# rem
# rem resamp() returns a filtered new sample point
#
# sub resamp(x, indat, alim, fmax, fsr, wnwdth)
# local i,j, r_g,r_w,r_a,r_snc,r_y : rem some local variables
# r_g = 2 * fmax / fsr : rem Calc gain correction factor
# r_y = 0
# for i = -wnwdth/2 to (wnwdth/2)-1 : rem For 1 window width
# j = int(x + i) : rem Calc input sample index
# : rem calculate von Hann Window. Scale and calculate Sinc
# r_w = 0.5 - 0.5 * cos(2*pi*(0.5 + (j - x)/wnwdth))
# r_a = 2*pi*(j - x)*fmax/fsr
# r_snc = 1 : if (r_a <> 0) then r_snc = sin(r_a)/r_a
# if (j >= 0) and (j < alim) then
# r_y = r_y + r_g * r_w * r_snc * indat(j)
# endif
# next i
# resamp = r_y : rem Return new filtered sample
# end sub
#
# rem - Ron Nicholson's QDSS ReSampler cookbook recipe
# rem QDSS = Quick, Dirty, Simple and Short
# rem Version 0.1b - 2007-Aug-01
# rem Copyright 2007 Ronald H. Nicholson Jr.
# rem No warranties implied. Error checking, optimization, and
# rem quality assessment of the "results" is left as an exercise
# rem for the student.
# rem (consider this code Open Source under a BSD style license)
# rem IMHO. YMMV. http://www.nicholson.com/rhn/dsp.html