Improvement to the Python wrapper
Corrected issues with arm_fir_decimate and arm_fir_interpolate Corrected issues with real FFTs in the wrapper Added a customization option for the FIFO class in compute graph. Added Python tests for the corrected functions.pull/67/head
Christophe Favergeon
3 years ago
parent
0705c67568
commit
3d1a7f7ff4
@ -0,0 +1,163 @@
|
||||
# Bug corrections for version 1.9
|
||||
import cmsisdsp as dsp
|
||||
import cmsisdsp.fixedpoint as f
|
||||
import numpy as np
|
||||
import math
|
||||
import colorama
|
||||
from colorama import init,Fore, Back, Style
|
||||
from numpy.testing import assert_allclose
|
||||
import matplotlib.pyplot as plt
|
||||
from scipy import signal
|
||||
|
||||
init()
|
||||
|
||||
def printTitle(s):
|
||||
print("\n" + Fore.GREEN + Style.BRIGHT + s + Style.RESET_ALL)
|
||||
|
||||
def printSubTitle(s):
|
||||
print("\n" + Style.BRIGHT + s + Style.RESET_ALL)
|
||||
|
||||
printTitle("Decimate")
|
||||
|
||||
test_length_seconds = 0.1
|
||||
signal_frequency = 100
|
||||
sampling_freq = 8000
|
||||
nbSamples = int(test_length_seconds*sampling_freq)
|
||||
wave = np.sin(2*np.pi*signal_frequency*np.linspace(0,test_length_seconds,nbSamples))
|
||||
|
||||
#plt.plot(wave)
|
||||
#plt.show()
|
||||
|
||||
decimationFactor = 4
|
||||
|
||||
numTaps = 9
|
||||
downsamplingFilter = signal.firwin(numTaps,1.0 / decimationFactor)
|
||||
block_size = 160
|
||||
assert(block_size % decimationFactor == 0)
|
||||
|
||||
ds_state = np.zeros(block_size + len(downsamplingFilter)-1)
|
||||
decimator = dsp.arm_fir_decimate_instance_f32()
|
||||
|
||||
status = dsp.arm_fir_decimate_init_f32(decimator,numTaps,decimationFactor, downsamplingFilter, ds_state)
|
||||
|
||||
def processSignal(sig,dec,f):
|
||||
result = []
|
||||
for blockNb in range(len(sig) // block_size):
|
||||
s = blockNb * block_size
|
||||
e = s + block_size
|
||||
|
||||
r = f(dec,sig[s:e])
|
||||
result.append(r)
|
||||
output = np.hstack(result)
|
||||
return(output)
|
||||
|
||||
ref = processSignal(wave,decimator,dsp.arm_fir_decimate_f32)
|
||||
#plt.plot(ref)
|
||||
#plt.show()
|
||||
|
||||
printSubTitle("Decimate Q31")
|
||||
waveQ31 = f.toQ31(wave)
|
||||
downsamplingFilterQ31 = f.toQ31(downsamplingFilter)
|
||||
stateQ31 = np.zeros(block_size + len(downsamplingFilter)-1)
|
||||
decimatorQ31 = dsp.arm_fir_decimate_instance_q31()
|
||||
status = dsp.arm_fir_decimate_init_q31(decimatorQ31,numTaps,decimationFactor,
|
||||
downsamplingFilterQ31, stateQ31)
|
||||
|
||||
outputQ31 = processSignal(waveQ31,decimatorQ31,dsp.arm_fir_decimate_q31)
|
||||
outputF32 = f.Q31toF32(outputQ31)
|
||||
|
||||
printSubTitle("Decimate Fast Q31")
|
||||
waveQ31 = f.toQ31(wave)
|
||||
downsamplingFilterQ31 = f.toQ31(downsamplingFilter)
|
||||
stateQ31 = np.zeros(block_size + len(downsamplingFilter)-1)
|
||||
decimatorQ31 = dsp.arm_fir_decimate_instance_q31()
|
||||
status = dsp.arm_fir_decimate_init_q31(decimatorQ31,numTaps,decimationFactor,
|
||||
downsamplingFilterQ31, stateQ31)
|
||||
|
||||
outputQ31 = processSignal(waveQ31,decimatorQ31,dsp.arm_fir_decimate_fast_q31)
|
||||
outputF32 = f.Q31toF32(outputQ31)
|
||||
|
||||
printSubTitle("Decimate Q15")
|
||||
|
||||
waveQ15 = f.toQ15(wave)
|
||||
downsamplingFilterQ15 = f.toQ15(downsamplingFilter)
|
||||
stateQ15 = np.zeros(block_size + len(downsamplingFilter)-1)
|
||||
decimatorQ15 = dsp.arm_fir_decimate_instance_q15()
|
||||
status = dsp.arm_fir_decimate_init_q15(decimatorQ15,numTaps,decimationFactor,
|
||||
downsamplingFilterQ15, stateQ15)
|
||||
|
||||
outputQ15 = processSignal(waveQ15,decimatorQ15,dsp.arm_fir_decimate_q15)
|
||||
outputF32 = f.Q15toF32(outputQ15)
|
||||
#plt.plot(outputF32)
|
||||
#plt.show()
|
||||
assert_allclose(ref,outputF32,rtol=2e-3,atol=1e-3)
|
||||
|
||||
printSubTitle("Decimate Fast Q15")
|
||||
|
||||
waveQ15 = f.toQ15(wave)
|
||||
downsamplingFilterQ15 = f.toQ15(downsamplingFilter)
|
||||
stateQ15 = np.zeros(block_size + len(downsamplingFilter)-1)
|
||||
decimatorQ15 = dsp.arm_fir_decimate_instance_q15()
|
||||
status = dsp.arm_fir_decimate_init_q15(decimatorQ15,numTaps,decimationFactor,
|
||||
downsamplingFilterQ15, stateQ15)
|
||||
|
||||
outputQ15 = processSignal(waveQ15,decimatorQ15,dsp.arm_fir_decimate_fast_q15)
|
||||
outputF32 = f.Q15toF32(outputQ15)
|
||||
#plt.plot(outputF32)
|
||||
#plt.show()
|
||||
assert_allclose(ref,outputF32,rtol=2e-3,atol=1e-3)
|
||||
|
||||
|
||||
printTitle("Interpolate")
|
||||
|
||||
upsamplingFactor = 4
|
||||
|
||||
numTaps = 16
|
||||
upsamplingFilter = signal.firwin(numTaps,1.0 / upsamplingFactor)
|
||||
assert(numTaps % upsamplingFactor == 0)
|
||||
block_size = 40
|
||||
|
||||
|
||||
printSubTitle("Interpolate F32")
|
||||
|
||||
state = np.zeros(block_size + len(upsamplingFilter)//upsamplingFactor-1)
|
||||
interpolator = dsp.arm_fir_interpolate_instance_f32()
|
||||
status = dsp.arm_fir_interpolate_init_f32(interpolator,upsamplingFactor,numTaps,
|
||||
upsamplingFilter, state)
|
||||
|
||||
output = processSignal(ref,interpolator,dsp.arm_fir_interpolate_f32)
|
||||
output = output / np.max(output)
|
||||
#t = range(nbSamples)
|
||||
#plt.plot(t,wave,t[:-11],output[11:])
|
||||
#plt.show()
|
||||
|
||||
d = 11
|
||||
assert_allclose(wave[:-d],output[d:],atol=0.1)
|
||||
|
||||
printSubTitle("Interpolate Q31")
|
||||
|
||||
upsamplingFilterQ31 = f.toQ31(upsamplingFilter)
|
||||
stateQ31 = np.zeros(block_size + len(upsamplingFilter)//upsamplingFactor-1)
|
||||
interpolatorQ31 = dsp.arm_fir_interpolate_instance_q31()
|
||||
status = dsp.arm_fir_interpolate_init_q31(interpolatorQ31,upsamplingFactor,numTaps,
|
||||
upsamplingFilterQ31, stateQ31)
|
||||
|
||||
outputQ31 = processSignal(outputQ31,interpolatorQ31,dsp.arm_fir_interpolate_q31)
|
||||
outputF32 = f.Q31toF32(outputQ31)
|
||||
outputF32 = outputF32 / np.max(outputF32)
|
||||
|
||||
assert_allclose(wave[:-d],outputF32[d:],atol=0.1)
|
||||
|
||||
printSubTitle("Interpolate Q15")
|
||||
|
||||
upsamplingFilterQ15 = f.toQ15(upsamplingFilter)
|
||||
stateQ15 = np.zeros(block_size + len(upsamplingFilter)//upsamplingFactor-1)
|
||||
interpolatorQ15 = dsp.arm_fir_interpolate_instance_q15()
|
||||
status = dsp.arm_fir_interpolate_init_q15(interpolatorQ15,upsamplingFactor,numTaps,
|
||||
upsamplingFilterQ15, stateQ15)
|
||||
|
||||
outputQ15 = processSignal(outputQ15,interpolatorQ15,dsp.arm_fir_interpolate_q15)
|
||||
outputF32 = f.Q15toF32(outputQ15)
|
||||
outputF32 = outputF32 / np.max(outputF32)
|
||||
|
||||
assert_allclose(wave[:-d],outputF32[d:],atol=0.1)
|
||||
@ -0,0 +1,146 @@
|
||||
import cmsisdsp as dsp
|
||||
import cmsisdsp.fixedpoint as f
|
||||
|
||||
import numpy as np
|
||||
from scipy import signal
|
||||
import matplotlib.pyplot as plt
|
||||
import scipy.fft
|
||||
|
||||
import colorama
|
||||
from colorama import init,Fore, Back, Style
|
||||
from numpy.testing import assert_allclose
|
||||
|
||||
init()
|
||||
|
||||
def printTitle(s):
|
||||
print("\n" + Fore.GREEN + Style.BRIGHT + s + Style.RESET_ALL)
|
||||
|
||||
def printSubTitle(s):
|
||||
print("\n" + Style.BRIGHT + s + Style.RESET_ALL)
|
||||
|
||||
|
||||
def chop(A, eps = 1e-6):
|
||||
B = np.copy(A)
|
||||
B[np.abs(A) < eps] = 0
|
||||
return B
|
||||
|
||||
nb = 32
|
||||
signal = np.cos(2 * np.pi * np.arange(nb) / nb)*np.cos(0.2*2 * np.pi * np.arange(nb) / nb)
|
||||
|
||||
ref=scipy.fft.rfft(signal)
|
||||
invref = scipy.fft.irfft(ref)
|
||||
|
||||
# Convert ref to CMSIS-DSP format
|
||||
referenceFloat=np.zeros(nb)
|
||||
# Replace complex datatype by real datatype
|
||||
referenceFloat[0::2] = np.real(ref)[:-1]
|
||||
referenceFloat[1::2] = np.imag(ref)[:-1]
|
||||
# Copy Nyquist frequency value into first
|
||||
# sample.This is just a storage trick so that the
|
||||
# output of the RFFT has same length as input
|
||||
# It is legacy behavior that we need to keep
|
||||
# for backward compatibility but it is not
|
||||
# very pretty
|
||||
referenceFloat[1] = np.real(ref[-1])
|
||||
|
||||
printTitle("RFFT FAST F64")
|
||||
|
||||
printSubTitle("RFFT")
|
||||
|
||||
|
||||
rfftf64=dsp.arm_rfft_fast_instance_f64()
|
||||
status=dsp.arm_rfft_fast_init_f64(rfftf64,nb)
|
||||
result = dsp.arm_rfft_fast_f64(rfftf64,signal,0)
|
||||
|
||||
|
||||
assert_allclose(referenceFloat,result)
|
||||
|
||||
printSubTitle("RIFFT")
|
||||
|
||||
rifftf64=dsp.arm_rfft_fast_instance_f64()
|
||||
status=dsp.arm_rfft_fast_init_f64(rifftf64,nb)
|
||||
result = dsp.arm_rfft_fast_f64(rifftf64,referenceFloat,1)
|
||||
|
||||
assert_allclose(invref,result,atol=1e-15)
|
||||
|
||||
printTitle("RFFT FAST F32")
|
||||
|
||||
printSubTitle("RFFT")
|
||||
|
||||
|
||||
rfftf32=dsp.arm_rfft_fast_instance_f32()
|
||||
status=dsp.arm_rfft_fast_init_f32(rfftf32,nb)
|
||||
result = dsp.arm_rfft_fast_f32(rfftf32,signal,0)
|
||||
|
||||
|
||||
assert_allclose(referenceFloat,result,rtol=3e-6)
|
||||
|
||||
printSubTitle("RIFFT")
|
||||
|
||||
rifftf32=dsp.arm_rfft_fast_instance_f32()
|
||||
status=dsp.arm_rfft_fast_init_f32(rifftf32,nb)
|
||||
result = dsp.arm_rfft_fast_f32(rifftf32,referenceFloat,1)
|
||||
|
||||
assert_allclose(invref,result,atol=1e-7)
|
||||
|
||||
# Fixed point
|
||||
# Reference from fixed point arithmetric.
|
||||
# The RFFT are not packing the Nyquist frequency
|
||||
# real value in sample 0
|
||||
referenceFloat=np.zeros(nb+2)
|
||||
# Replace complex datatype by real datatype
|
||||
referenceFloat[0::2] = np.real(ref)
|
||||
referenceFloat[1::2] = np.imag(ref)
|
||||
|
||||
printTitle("RFFT Q31")
|
||||
|
||||
printSubTitle("RFFT")
|
||||
|
||||
signalQ31 = f.toQ31(signal)
|
||||
rfftQ31=dsp.arm_rfft_instance_q31()
|
||||
status=dsp.arm_rfft_init_q31(rfftQ31,nb,0,1)
|
||||
resultQ31 = dsp.arm_rfft_q31(rfftQ31,signalQ31)
|
||||
# Drop the conjugate part which is not computed by scipy
|
||||
resultQ31 = resultQ31[:nb+2]
|
||||
resultF = f.Q31toF32(resultQ31) * nb
|
||||
|
||||
assert_allclose(referenceFloat,resultF,rtol=1e-6,atol=1e-6)
|
||||
|
||||
|
||||
printSubTitle("RIFFT")
|
||||
|
||||
rifftQ31=dsp.arm_rfft_instance_q31()
|
||||
status=dsp.arm_rfft_init_q31(rifftQ31,nb,1,1)
|
||||
# Apply CMSIS-DSP scaling
|
||||
referenceQ31 = f.toQ31(referenceFloat / nb)
|
||||
resultQ31 = dsp.arm_rfft_q31(rifftQ31,referenceFloat)
|
||||
resultF = f.Q31toF32(resultQ31)
|
||||
|
||||
assert_allclose(invref,result,atol=1e-6)
|
||||
|
||||
printTitle("RFFT Q15")
|
||||
|
||||
printSubTitle("RFFT")
|
||||
|
||||
signalQ15 = f.toQ15(signal)
|
||||
rfftQ15=dsp.arm_rfft_instance_q15()
|
||||
status=dsp.arm_rfft_init_q15(rfftQ15,nb,0,1)
|
||||
resultQ15 = dsp.arm_rfft_q15(rfftQ15,signalQ15)
|
||||
# Drop the conjugate part which is not computed by scipy
|
||||
resultQ15 = resultQ15[:nb+2]
|
||||
resultF = f.Q15toF32(resultQ15) * nb
|
||||
|
||||
assert_allclose(referenceFloat,resultF,rtol=1e-6,atol=1e-2)
|
||||
|
||||
|
||||
printSubTitle("RIFFT")
|
||||
|
||||
rifftQ15=dsp.arm_rfft_instance_q15()
|
||||
status=dsp.arm_rfft_init_q15(rifftQ15,nb,1,1)
|
||||
# Apply CMSIS-DSP scaling
|
||||
referenceQ15 = f.toQ15(referenceFloat / nb)
|
||||
resultQ15 = dsp.arm_rfft_q15(rifftQ15,referenceFloat)
|
||||
resultF = f.Q15toF32(resultQ15)
|
||||
|
||||
assert_allclose(invref,result,atol=1e-2)
|
||||
|
||||
@ -1,27 +0,0 @@
|
||||
import cmsisdsp as dsp
|
||||
import numpy as np
|
||||
from scipy import signal
|
||||
import matplotlib.pyplot as plt
|
||||
import scipy.fft
|
||||
|
||||
|
||||
def chop(A, eps = 1e-6):
|
||||
B = np.copy(A)
|
||||
B[np.abs(A) < eps] = 0
|
||||
return B
|
||||
|
||||
nb = 32
|
||||
signal = np.cos(2 * np.pi * np.arange(nb) / nb)*np.cos(0.2*2 * np.pi * np.arange(nb) / nb)
|
||||
|
||||
#print("{")
|
||||
#for x in signal:
|
||||
# print("%f," % x)
|
||||
#print("}")
|
||||
|
||||
result1=scipy.fft.rfft(signal)
|
||||
print(chop(result1))
|
||||
rfftf32=dsp.arm_rfft_fast_instance_f32()
|
||||
status=dsp.arm_rfft_fast_init_f32(rfftf32,nb)
|
||||
print(status)
|
||||
resultI = dsp.arm_rfft_fast_f32(rfftf32,signal,0)
|
||||
print(chop(resultI))
|
||||
@ -1,2 +1,2 @@
|
||||
# Python wrapper version
|
||||
__version__ = "1.9.2"
|
||||
__version__ = "1.9.3"
|
||||
|
||||
Loading…
Reference in New Issue