Permalink
Browse files

KFR 3.0.2

  • Loading branch information...
dlevin256 committed Dec 19, 2018
1 parent 39488b8 commit b73a4f29e22f935e53accfd5ccefe276af0bf0de
@@ -1,3 +1,6 @@
# Backup
*.bak

# Compiled Object files
*.slo
*.lo
@@ -29,7 +29,7 @@ def gen_tick_labels(stop, start=10):
yield ''
for t in gen_tick_labels(stop, start * 10):
yield t

def smooth_colormap(colors, name='cmap1'):
to_rgb = clr.ColorConverter().to_rgb
colors = [(p, to_rgb(c)) for p, c in colors]
@@ -57,21 +57,21 @@ def wavplot(data, title='Title', file=None, segmentsize=512, overlap=8, Fs=48000
(7/9, '#fdc967'),
(8/9, '#f3fab8'),
(1 , '#ffffff')
])
])

if not isinstance(data, (list, tuple, np.ndarray)):
w = wave.open(data, 'rb')
Fs = w.getframerate()
data = np.fromstring(w.readframes(w.getnframes()), dtype=np.int32)/2147483647.0
data = np.fromstring(w.readframes(w.getnframes()), dtype=np.int32)/2147483647.0
data = np.array(data)
if normalize:
maxitem = max(abs(np.max(data)), abs(np.min(data)))
if maxitem > 0:
data = data / maxitem

datalen = len(data)

def fast_resample(data, newlen):
def fast_resample(data, newlen):
oldlen=len(data)
result=[]
for i in range(newlen):
@@ -92,24 +92,24 @@ def fast_resample(data, newlen):
r = range(segm*datalen//segments, segm*datalen//segments+segmentsize*overlap)
subdata = data[r]
subdata = subdata * window
n = len(subdata)
n = len(subdata)
Y = np.fft.fft(subdata)/n
Y = Y[range(len(Y) // 2)]
Yfreq = 20 * np.log10(np.absolute(Y))
Yfreq = signal.resample(Yfreq, 512)
Yfreq = np.fmax(-300, Yfreq)
im.append(Yfreq)

im = np.transpose(im)

plt.imshow(im,cmap=cmap, aspect='auto', vmin=vmin, vmax=vmax, origin='lower', extent=[0, datalen / Fs, 0, Fs / 2 ], interpolation='bicubic')
plt.colorbar()

if not file:
plt.show()
else:
plt.savefig(file)


def plot(data,
title='Title',
@@ -128,7 +128,15 @@ def plot(data,
spectrogram=False,
vmin=-160,
vmax=0,
normalize=False):
normalize=False,
freqticks=[],
freq_lim=None,
freq_dB_lim=None,
phasearg=None):

if phasearg == 'auto':
phasearg = (len(data) - 1) * 0.5

if isinstance(data, (list, tuple, np.ndarray)) and not spectrogram:
if normalize:
maxitem = max(abs(np.max(data)), abs(np.min(data)))
@@ -154,7 +162,7 @@ def plot(data,
dataplot.set_autoscalex_on(False)
dataplot.set_xlim([0, n - 1])
dataplot.set_ylim(bottom=np.min(data))

np.seterr(all='ignore')

if freqresp or phaseresp:
@@ -180,11 +188,17 @@ def set_freq(a):
a.set_xticks(list(gen_ticks(Fs / 2)))
a.set_xticklabels(list(gen_tick_labels(Fs / 2)))
X = np.linspace(0, Fs / 2, len(Y), False)
a.set_xlim([10, Fs / 2])
if freq_lim is not None:
a.set_xlim(freq_lim)
else:
a.set_xlim([10, Fs / 2])
a.set_xticks(list(a.get_xticks()) + freqticks)
return X

if freqresp:
freqplot = a[1]
if freq_dB_lim is not None:
freqplot.set_ylim(freq_dB_lim)
X = set_freq(freqplot)
freqplot.set_ylabel('Gain (dB)')
freqplot.grid(True, **grid_style)
@@ -193,12 +207,15 @@ def set_freq(a):

if phaseresp:
phaseplot = a[1 + freqresp]
Yphase = np.angle(Y, deg=True);
if phasearg is not None:
Y = Y / 1j ** np.linspace(0, -(phasearg * 2), len(Y), endpoint=False)
Yphase = np.angle(Y, deg=True)
Yphase = np.select([Yphase < -179, True], [Yphase + 360, Yphase])
X = set_freq(phaseplot)
phaseplot.grid(True, **grid_style)
phaseplot.set_ylabel(r'Phase (${\circ}$)')
phaseplot.set_ylabel(r'Phase (${\circ}$)' if phasearg is None else r'Phase shift (${\circ}$)')
phaseplot.set_autoscaley_on(False)
phaseplot.set_ylim([-180, +180])
phaseplot.set_ylim([-190, +190])
phaseplot.plot(X, Yphase, **style)

plt.tight_layout(rect=[0, 0.0, 1, 0.94])
@@ -210,15 +227,15 @@ def set_freq(a):
else:
wavplot(data, title=title, file=file, segmentsize=segmentsize, overlap=overlap, Fs=Fs, vmin=vmin, vmax=vmax, normalize=normalize)


def perfplot(data, labels, title='Speed', xlabel='X', units='ms', file=None):

styles = [
{'color': '#F6511D', 'linestyle': '-', 'marker': 'o', 'markersize': 10.0, 'markeredgecolor': '#FFFFFF'},
{'color': '#00A6ED', 'linestyle': '-', 'marker': 'o', 'markersize': 10.0, 'markeredgecolor': '#FFFFFF'},
{'color': '#FFB400', 'linestyle': '-', 'marker': 'o', 'markersize': 10.0, 'markeredgecolor': '#FFFFFF'},
{'color': '#7FB800', 'linestyle': '-', 'marker': 'o', 'markersize': 10.0, 'markeredgecolor': '#FFFFFF'},
{'color': '#0D2C54', 'linestyle': '-', 'marker': 'o', 'markersize': 10.0, 'markeredgecolor': '#FFFFFF'},
{'color': '#0D2C54', 'linestyle': '-', 'marker': 'o', 'markersize': 10.0, 'markeredgecolor': '#FFFFFF'},
]
grid_style = {'color': '#777777'}
fig, ax = plt.subplots()
@@ -230,14 +247,14 @@ def perfplot(data, labels, title='Speed', xlabel='X', units='ms', file=None):
ax.set_xlabel(xlabel)
ax.set_ylabel(units)
x = np.linspace(0,len(d),len(d), False)
ax.plot(x, d, linewidth=1.6, label=l, **s)
ax.plot(x, d, linewidth=1.6, label=l, **s)

ax.set_ylim(bottom=0.0)
legend = ax.legend(loc='lower center', shadow=True)

plt.xticks(x, ticks, rotation='vertical')
plt.tight_layout(rect=[0, 0.0, 1, 0.94])

if not file:
plt.show()
else:
@@ -28,7 +28,7 @@ add_executable(window window.cpp)
target_link_libraries(window kfr)

add_executable(fir fir.cpp)
target_link_libraries(fir kfr)
target_link_libraries(fir kfr kfr_dft)

add_executable(sample_rate_conversion sample_rate_conversion.cpp)
target_link_libraries(sample_rate_conversion kfr kfr_io)
@@ -5,49 +5,144 @@
*/

#include <kfr/base.hpp>
#include <kfr/dft.hpp>
#include <kfr/dsp.hpp>
#include <kfr/io.hpp>

using namespace kfr;

#ifndef PYTHON_IS_INSTALLED
#define PYTHON_IS_INSTALLED 1
#endif

int main()
{
// Show KFR version
println(library_version());

const std::string options = "phaseresp=False";
// --------------------------------------------------------------------------------------
// --------------------- FIR filter design (using window functions) ---------------------
// --------------------------------------------------------------------------------------

// Options for dspplot
const std::string options = "phaseresp=True";

univector<fbase, 15> taps15;
univector<fbase, 127> taps127;
univector<fbase, 8191> taps8191;

// Prepare window functions (only expression saved here, not data)
expression_pointer<fbase> hann = to_pointer(window_hann(taps15.size()));

expression_pointer<fbase> kaiser = to_pointer(window_kaiser(taps127.size(), 3.0));

expression_pointer<fbase> blackman_harris = to_pointer(window_blackman_harris(taps8191.size()));

// Fill taps15 with the low pass FIR filter coefficients using hann window and cutoff=0.15
fir_lowpass(taps15, 0.15, hann, true);
plot_save("fir_lowpass_hann", taps15, options + ", title='15-point lowpass FIR, Hann window'");
#if PYTHON_IS_INSTALLED
// Plot filter, frequency and impulse response
// plot_save calls python (matplotlib and numpy must be installed) and saves SVG file
plot_save("fir_lowpass_hann", taps15,
options + ", phasearg='auto', title='15-point lowpass FIR, Hann window'");
#endif

// Fill taps127 with the low pass FIR filter coefficients using kaiser window and cutoff=0.2
fir_lowpass(taps127, 0.2, kaiser, true);
#if PYTHON_IS_INSTALLED
// Plot filter, frequency and impulse response
plot_save("fir_lowpass_kaiser", taps127,
options + ", title=r'127-point lowpass FIR, Kaiser window ($\\alpha=3.0$)'");
options + ", phasearg='auto', title=r'127-point lowpass FIR, Kaiser window ($\\alpha=3.0$)'");
#endif

// Fill taps127 with the high pass FIR filter coefficients using kaiser window and cutoff=0.2
fir_highpass(taps127, 0.2, kaiser, true);
#if PYTHON_IS_INSTALLED
// Plot filter, frequency and impulse response
plot_save("fir_highpass_kaiser", taps127,
options + ", title=r'127-point highpass FIR, Kaiser window ($\\alpha=3.0$)'");
options + ", phasearg='auto', title=r'127-point highpass FIR, Kaiser window ($\\alpha=3.0$)'");
#endif

// Fill taps127 with the band pass FIR filter coefficients using kaiser window and cutoff=0.2 and 0.4
fir_bandpass(taps127, 0.2, 0.4, kaiser, true);
#if PYTHON_IS_INSTALLED
// Plot filter, frequency and impulse response
plot_save("fir_bandpass_kaiser", taps127,
options + ", title=r'127-point bandpass FIR, Kaiser window ($\\alpha=3.0$)'");
options + ", phasearg='auto', title=r'127-point bandpass FIR, Kaiser window ($\\alpha=3.0$)'");
#endif

// Fill taps127 with the band stop FIR filter coefficients using kaiser window and cutoff=0.2 and 0.4
fir_bandstop(taps127, 0.2, 0.4, kaiser, true);
#if PYTHON_IS_INSTALLED
// Show filter, frequency and impulse response
plot_save("fir_bandstop_kaiser", taps127,
options + ", title=r'127-point bandstop FIR, Kaiser window ($\\alpha=3.0$)'");
options + ", phasearg='auto', title=r'127-point bandstop FIR, Kaiser window ($\\alpha=3.0$)'");
#endif

// Fill taps8191 with the low pass FIR filter coefficients using blackman harris window and cutoff=0.15
fir_lowpass(taps8191, 0.15, blackman_harris, true);
plot_save("fir_lowpass_blackman", taps8191,
options + ", title='8191-point lowpass FIR, Blackman-Harris window'");
univector<fbase, 8191 + 150> taps8191_150 = scalar(0);

// Shift by 150 samples
taps8191_150.slice(150) = taps8191;

#if PYTHON_IS_INSTALLED
// Plot filter, frequency and impulse response, pass phasearg to get correct phase shift (phasearg=offset
// to unit impulse in samples)
plot_save(
"fir_lowpass_blackman", taps8191_150,
options +
", phasearg=4095+150, title='8191-point lowpass FIR, Blackman-Harris window', padwidth=16384");
#endif

// --------------------------------------------------------------------------------------
// -------------------------- Using FIR filter as an expression -------------------------
// --------------------------------------------------------------------------------------

// Prepare 10000 samples of white noise
univector<float> noise = truncate(gen_random_range(-1.f, +1.f), 10000);

// Apply band stop filter
univector<float> filtered_noise = fir(noise, taps127);

#if PYTHON_IS_INSTALLED
// Plot results
plot_save("noise", noise, "title='Original noise', div_by_N=True");
plot_save("filtered_noise", filtered_noise, "title='Filtered noise', div_by_N=True");
#endif

// --------------------------------------------------------------------------------------
// ------------------------------- FIR filter as a class --------------------------------
// --------------------------------------------------------------------------------------

fir_bandpass(taps127, 0.2, 0.4, kaiser, true);
// Initialize FIR filter with float input/output and fbase taps
filter_fir<fbase, float> fir_filter(taps127);

// Apply to univector, static array, data by pointer or anything
univector<float> filtered_noise2;
fir_filter.apply(filtered_noise2, noise);

#if PYTHON_IS_INSTALLED
// Plot results
plot_save("filtered_noise2", filtered_noise2, "title='Filtered noise 2', div_by_N=True");
#endif

// --------------------------------------------------------------------------------------
// ---------------------- Convolution filter (optimized using DFT) ----------------------
// --------------------------------------------------------------------------------------

// Initialize FIR filter with float input/output and fbase taps
convolve_filter<float> conv_filter(taps127);

// Apply to univector, static array, data by pointer or anything
univector<float> filtered_noise3;
conv_filter.apply(filtered_noise3, noise);

#if PYTHON_IS_INSTALLED
// Plot results, same as filtered_noise2
plot_save("filtered_noise3", filtered_noise3, "title='Filtered noise 3', div_by_N=True");
#endif

return 0;
}
@@ -18,12 +18,10 @@ int main()
{
println(library_version());

const std::string options = "phaseresp=False";

univector<fbase> swept_sine = swept(0.5, len);

{
auto r = resampler<fbase>(resample_quality::high, output_sr, input_sr, 1.0, 0.496);
auto r = resampler<fbase>(resample_quality::high, output_sr, input_sr);
univector<fbase> resampled(len * output_sr / input_sr + r.get_delay());
r.process(resampled, swept_sine);

@@ -36,7 +34,7 @@ int main()
}

{
auto r = resampler<fbase>(resample_quality::normal, output_sr, input_sr, 1.0, 0.496);
auto r = resampler<fbase>(resample_quality::normal, output_sr, input_sr);
univector<fbase> resampled(len * output_sr / input_sr + r.get_delay());
r.process(resampled, swept_sine);

@@ -49,7 +47,7 @@ int main()
}

{
auto r = resampler<fbase>(resample_quality::low, output_sr, input_sr, 1.0, 0.496);
auto r = resampler<fbase>(resample_quality::low, output_sr, input_sr);
univector<fbase> resampled(len * output_sr / input_sr + r.get_delay());
r.process(resampled, swept_sine);

@@ -62,7 +60,7 @@ int main()
}

{
auto r = resampler<fbase>(resample_quality::draft, output_sr, input_sr, 1.0, 0.496);
auto r = resampler<fbase>(resample_quality::draft, output_sr, input_sr);
univector<fbase> resampled(len * output_sr / input_sr + r.get_delay());
r.process(resampled, swept_sine);

Oops, something went wrong.

0 comments on commit b73a4f2

Please sign in to comment.