219 lines
7.5 KiB
Python
219 lines
7.5 KiB
Python
import unittest
|
|
import numpy as np
|
|
from tests.utils_testing import get_path_for_data_file
|
|
from urh.signalprocessing.Filter import Filter
|
|
from urh.signalprocessing.Modulator import Modulator
|
|
from urh.signalprocessing.Signal import Signal
|
|
import array
|
|
|
|
from matplotlib import pyplot as plt
|
|
from urh.cythonext import signal_functions
|
|
from urh.signalprocessing.Spectrogram import Spectrogram
|
|
|
|
|
|
class SpectrogramTest(unittest.TestCase):
|
|
""" short time fourier transform of audio signal """
|
|
|
|
def stft(self, samples, window_size, overlap_factor=0.5, window_function=np.hanning):
|
|
"""
|
|
Perform Short-time Fourier transform to get the spectrogram for the given samples
|
|
|
|
:param samples: Complex samples
|
|
:param window_size: Size of DFT window
|
|
:param overlap_factor: Value between 0 (= No Overlapping) and 1 (= Full overlapping) of windows
|
|
:param window_function: Function for DFT window
|
|
:return: short-time Fourier transform of the given signal
|
|
"""
|
|
window = window_function(window_size)
|
|
|
|
# hop size determines by how many samples the window is advanced
|
|
hop_size = window_size - int(overlap_factor * window_size)
|
|
|
|
# pad with zeros to ensure last window fits signal
|
|
padded_samples = np.append(samples, np.zeros((len(samples) - window_size) % hop_size))
|
|
num_frames = ((len(padded_samples) - window_size) // hop_size) + 1
|
|
frames = [padded_samples[i*hop_size:i*hop_size+window_size] * window for i in range(num_frames)]
|
|
return np.fft.fft(frames)
|
|
|
|
def setUp(self):
|
|
self.signal = Signal(get_path_for_data_file("two_participants.complex16s"), "test")
|
|
|
|
def test_numpy_impl(self):
|
|
sample_rate = 1e6
|
|
spectrogram = np.fft.fftshift(self.stft(self.signal.iq_array.data, 2**10, overlap_factor=0.5)) / 1024
|
|
|
|
ims = 10 * np.log10(spectrogram.real ** 2 + spectrogram.imag ** 2) # convert amplitudes to decibel
|
|
num_time_bins, num_freq_bins = np.shape(ims)
|
|
|
|
plt.imshow(np.transpose(ims), aspect="auto", cmap="magma")
|
|
plt.colorbar()
|
|
|
|
plt.xlabel("time in seconds")
|
|
plt.ylabel("frequency in Hz")
|
|
plt.ylim(ymin=0, ymax=num_freq_bins)
|
|
|
|
x_tick_pos = np.linspace(0, num_time_bins - 1, 5, dtype=np.float32)
|
|
plt.xticks(x_tick_pos, ["%.02f" % l for l in (x_tick_pos * len(self.signal.iq_array.data) / num_time_bins) / sample_rate])
|
|
y_tick_pos = np.linspace(0, num_freq_bins - 1, 10, dtype=np.int16)
|
|
frequencies = np.fft.fftshift(np.fft.fftfreq(num_freq_bins, 1/sample_rate))
|
|
plt.yticks(y_tick_pos, ["%.02f" % frequencies[i] for i in y_tick_pos])
|
|
|
|
plt.show()
|
|
|
|
def narrowband_iir(self, fc, bw, fs):
|
|
fc /= fs
|
|
bw /= fs
|
|
|
|
R = 1 - 3 * bw
|
|
K = (1 - 2 * R * np.cos(2 * np.pi * fc) + R ** 2) / (2 - 2*np.cos(2 * np.pi * fc))
|
|
|
|
a = np.array([K, -2*K*np.cos(2 * np.pi * fc), K], dtype=np.float64)
|
|
b = np.array([2 * R * np.cos(2 * np.pi * fc), -R**2], dtype=np.float64)
|
|
|
|
return a, b
|
|
|
|
def test_bandpass(self):
|
|
# Generate a noisy signal
|
|
fs = 2000
|
|
T = 0.1
|
|
nsamples = T * fs
|
|
t = np.linspace(0, T, nsamples, endpoint=False)
|
|
a = 0.02
|
|
f0 = 600
|
|
x = 0.25 * np.sin(2 * np.pi * 0.25*f0 * t)
|
|
x += 0.25 * np.sin(2 * np.pi * f0 * t)
|
|
x += 0.25 * np.sin(2 * np.pi * 2*f0 * t)
|
|
x += 0.25 * np.sin(2 * np.pi * 3*f0 * t)
|
|
|
|
import time
|
|
|
|
lowcut = f0 - 200
|
|
highcut = f0 + 200
|
|
|
|
# Define the parameters
|
|
fc = f0 / fs
|
|
b = 0.05
|
|
data = x
|
|
|
|
y = Filter.apply_bandpass_filter(data, lowcut / fs, highcut / fs, filter_bw=b)
|
|
|
|
plt.plot(y, label='Filtered signal (%g Hz)' % f0)
|
|
plt.plot(data, label='Noisy signal')
|
|
plt.legend(loc='upper left')
|
|
plt.show()
|
|
|
|
def test_iir_bandpass(self):
|
|
# Generate a noisy signal
|
|
fs = 2400
|
|
T = 6
|
|
nsamples = T * fs
|
|
t = np.linspace(0, T, nsamples, endpoint=False)
|
|
a = 0.02
|
|
f0 = 300
|
|
x = 0.5 * np.sin(2 * np.pi * f0 * t)
|
|
x += 0.25 * np.sin(2 * np.pi * 2 * f0 * t)
|
|
x += 0.25 * np.sin(2 * np.pi * 3 * f0 * t)
|
|
|
|
#data = x.astype(np.complex64)
|
|
data = np.sin(2 * np.pi * f0 * t).astype(np.complex64)
|
|
|
|
print("Len data", len(data))
|
|
a, b = self.narrowband_iir(f0, 100, fs)
|
|
s = a.sum() + b.sum()
|
|
#a /= s
|
|
#b /= s
|
|
print(a, b)
|
|
|
|
filtered_data = signal_functions.iir_filter(a, b, data)
|
|
|
|
#plt.plot(data, label='Noisy signal')
|
|
plt.plot(np.fft.fft(filtered_data), label='Filtered signal (%g Hz)' % f0)
|
|
|
|
|
|
|
|
plt.legend(loc='upper left')
|
|
plt.show()
|
|
|
|
def test_channels(self):
|
|
sample_rate = 10 ** 6
|
|
|
|
channel1_freq = 40 * 10 ** 3
|
|
channel2_freq = 240 * 10 ** 3
|
|
|
|
channel1_data = array.array("B", [1, 0, 1, 0, 1, 0, 0, 1])
|
|
channel2_data = array.array("B", [1, 1, 0, 0, 1, 1, 0, 1])
|
|
channel3_data = array.array("B", [1, 0, 0, 1, 0, 1, 1, 1])
|
|
|
|
filter_bw = 0.1
|
|
|
|
filter_freq1_high = 1.5 * channel1_freq
|
|
filter_freq1_low = 0.5 * channel1_freq
|
|
filter_freq2_high = 1.5*channel2_freq
|
|
filter_freq2_low = 0.5 * channel2_freq
|
|
|
|
modulator1, modulator2, modulator3 = Modulator("test"), Modulator("test2"), Modulator("test3")
|
|
modulator1.carrier_freq_hz = channel1_freq
|
|
modulator2.carrier_freq_hz = channel2_freq
|
|
modulator3.carrier_freq_hz = -channel2_freq
|
|
modulator1.sample_rate = modulator2.sample_rate = modulator3.sample_rate = sample_rate
|
|
data1 = modulator1.modulate(channel1_data)
|
|
data2 = modulator2.modulate(channel2_data)
|
|
data3 = modulator3.modulate(channel3_data)
|
|
|
|
mixed_signal = data1 + data2 + data3
|
|
|
|
mixed_signal.tofile("/tmp/three_channels.complex")
|
|
|
|
plt.subplot("221")
|
|
plt.title("Signal")
|
|
plt.plot(mixed_signal)
|
|
|
|
spectrogram = Spectrogram(mixed_signal)
|
|
plt.subplot("222")
|
|
plt.title("Spectrogram")
|
|
plt.imshow(np.transpose(spectrogram.data), aspect="auto", cmap="magma")
|
|
plt.ylim(0, spectrogram.freq_bins)
|
|
|
|
chann1_filtered = Filter.apply_bandpass_filter(mixed_signal, filter_freq1_low / sample_rate, filter_freq1_high / sample_rate, filter_bw)
|
|
plt.subplot("223")
|
|
plt.title("Channel 1 Filtered ({})".format("".join(map(str, channel1_data))))
|
|
plt.plot(chann1_filtered)
|
|
|
|
chann2_filtered = Filter.apply_bandpass_filter(mixed_signal, filter_freq2_low / sample_rate, filter_freq2_high / sample_rate, filter_bw)
|
|
plt.subplot("224")
|
|
plt.title("Channel 2 Filtered ({})".format("".join(map(str, channel2_data))))
|
|
plt.plot(chann2_filtered)
|
|
|
|
plt.show()
|
|
|
|
|
|
def test_bandpass_h(self):
|
|
f_low = -0.4
|
|
f_high = -0.3
|
|
bw = 0.01
|
|
|
|
f_shift = (f_low + f_high) / 2
|
|
f_c = (f_high - f_low) / 2
|
|
|
|
N = Filter.get_filter_length_from_bandwidth(bw)
|
|
|
|
h = Filter.design_windowed_sinc_lpf(f_c, bw=bw) * np.exp(np.complex(0,1) * np.pi * 2 * f_shift * np.arange(0, N, dtype=complex))
|
|
|
|
#h = Filter.design_windowed_sinc_bandpass(f_low=f_low, f_high=f_high, bw=bw)
|
|
#h = Filter.design_windowed_sinc_lpf(0.42, bw=0.08)
|
|
|
|
impulse = np.exp(1j * np.linspace(0, 1, 50))
|
|
|
|
plt.subplot("221")
|
|
plt.title("f_low={} f_high={} bw={}".format(f_low, f_high, bw))
|
|
plt.plot(np.fft.fftfreq(1024), np.fft.fft(h, 1024))
|
|
|
|
plt.subplot("222")
|
|
plt.plot(h)
|
|
|
|
plt.show()
|
|
|
|
|
|
|
|
# h = cls.design_windowed_sinc_bandpass(f_low, f_high, filter_bw)
|