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)