Commits (4)
%% Cell type:code id: tags:
``` python
# Copyright 2019 Institut für Nachrichtentechnik, RWTH Aachen University
%matplotlib widget
from ient_plots import *
from ient_transforms import *
```
%% Cell type:markdown id: tags:
<div>
<img src="figures/rwth_ient_logo@2x.png" style="float: right;height: 5em;">
</div>
# Transforms
## Contents
Following transforms are defined in `ient_transforms`:
``` python
ient_dft(s, fs, NFFT)
ient_idft(S, Ntime, NFFT)
ient_ilaplace_ht(t, H0, pp, pz, ord_p, ord_z, roc)
ient_ilaplace_Hf(f, H0, pp, pz, ord_p, ord_z, dB)
ient_iz_ht(n, H0, pp, pz, ord_p, ord_z, roc)
ient_iz_Hf(f, H0, pp, pz, ord_p, ord_z, dB)
```
%% Cell type:markdown id: tags:
## Fourier Transform
%% Cell type:code id: tags:
``` python
# Time Domain
fs = 44100 # very high sampling rate assumed, to simulate quasi-continuous time and frequency axis
t = np.linspace(-2.5, 2.5, 5*fs)
s = np.sin(2*np.pi*500*t)
# Fourier Transform
S,f = ient_dft(s, fs)
fig,axs = plt.subplots(2,1, figsize=(8,6));
ax = axs[0]; ax.plot(t*1000, s, 'rwth');
ax.set_xlabel(r'$\rightarrow t$ [ms]', bbox=ient_wbbox); ax.set_ylabel(r'$\uparrow s(t)$', bbox=ient_wbbox)
ax.set_xlim([-11, 11]); ax.set_ylim([-1.1, 1.19]); ient_axis(ax);
ax = axs[1]; ax.plot(f, np.abs(S), 'rwth');
ax.set_xlabel(r'$\rightarrow f$ [Hz]', bbox=ient_wbbox); ax.set_ylabel(r'$\uparrow |S(f)|$', bbox=ient_wbbox)
ax.set_xlim([-1100, 1100]); ax.set_ylim([0, 0.65]); ient_axis(ax);
```
%% Cell type:markdown id: tags:
Inverse Fourier transform
%% Cell type:code id: tags:
``` python
s2 = ient_idft(S, len(s));
fig,ax = plt.subplots(1,1, figsize=(8,4));
ax.plot(t*1000, np.real(s2), 'rwth');
ax.set_xlabel(r'$\rightarrow t$ [ms]', bbox=ient_wbbox); ax.set_ylabel(r'$\uparrow \mathcal{F}^{-1}\{S(f)\}$', bbox=ient_wbbox)
ax.set_xlim([-11, 11]); ax.set_ylim([-1.1, 1.19]); ient_axis(ax);
```
%% Cell type:markdown id: tags:
## Laplace Transform
%% Cell type:markdown id: tags:
Pole-zero plot is explained in [Plots](Plots.ipynb)
%% Cell type:markdown id: tags:
Inverse Laplace Transform
%% Cell type:code id: tags:
``` python
fig,axs = plt.subplots(1,2, figsize=(10, 4))
t = np.linspace(-6, 6, 1024)
f = np.linspace(-6, 6, 1024)
pp = np.array([-2]); pz = np.array([]) # Poles and Zeros
ord_p = np.array([1]); ord_z = np.array([]) # Poles' and Zeros' orders
roc = np.array([-2, np.inf]) # region of convergence
H0 = 1
# Time Domain
s1, t1d , s1d = ient_ilaplace_ht(t, H0, pp, pz, ord_p, ord_z, roc)
axs[0].set_xlabel(r'$\rightarrow t$'); axs[0].set_ylabel(r'$\uparrow s_1(t)$'); ient_grid(axs[0]); ient_axis(axs[0])
axs[0].set_xlim([-5.5,5.5]); axs[0].set_ylim([-0.1,1.05]);
axs[0].plot(t, np.real(s1), **ient_style_graph)
ient_plot_dirac(axs[0], t1d, s1d);
# Frequency Domain
S1f = ient_ilaplace_Hf(f, H0, pp, pz, ord_p, ord_z, dB=False)
axs[1].set_xlabel(r'$\rightarrow f$'); axs[1].set_ylabel(r'$\uparrow S_1(f)$'); ient_grid(axs[1]); ient_axis(axs[1])
axs[1].set_xlim([-5.5,5.5]); axs[1].set_ylim([-0.1,0.55]);
axs[1].plot(f, S1f, **ient_style_graph);
```
%% Cell type:markdown id: tags:
This code is licensed under the [MIT license](https://opensource.org/licenses/MIT).
# z-Transform
%% Cell type:markdown id: tags:
Pole-zero plot is explained in [Plots](Plots.ipynb)
%% Cell type:code id: tags:
``` python
fig,axs = plt.subplots(1,2, figsize=(10, 4))
n = np.linspace(-6, 6, 13)
f = np.linspace(-6, 6, 1024)
H0 = 1
pp = np.array([.5, 2]); pz = np.array([0])
ord_p = np.array([1, 1]); ord_z = np.array([1])
roc = np.array([.5, 2])
# Discrete Time Domain
s_n = ient_iz_ht(n, H0, pp, pz, ord_p, ord_z, roc)
ax = axs[0]
ax.set_xlabel(r'$\rightarrow n$'); ax.set_ylabel(r'$\uparrow s(n)$'); ient_grid(ax); ient_axis(ax)
ax.set_xlim(np.min(n), np.max(n)); ient_update_ylim(ax, s_n, 0.19, 1e05);
ient_stem(ax, n, s_n)
# Frequency Domain
S_f = ient_iz_Hf(f, H0, pp, pz, ord_p, ord_z, dB=False)
ax = axs[1]
ax.set_xlabel(r'$\rightarrow f$'); ax.set_ylabel(r'$\uparrow S(f)$'); ient_grid(ax); ient_axis(ax)
ax.set_xlim([-5.5,5.5]); ient_update_ylim(ax, S_f, 0.19, 1e05);
ax.plot(f, S_f, **ient_style_graph);
```
%% Cell type:markdown id: tags:
This code is licensed under the [MIT license](https://opensource.org/licenses/MIT).
......
......@@ -15,16 +15,48 @@ rl_hp = lambda f: (1j * 2 * np.pi * f) / (1 + 1j * 2 * np.pi * f) # RL = 1
eps = np.finfo(float).eps # floating point accuracy
def findIndOfLeastDiff(x, x0):
def find_intervals(s, thresh, delta):
"""
Find Index of the value that's nearest to x0
Find intervals of signal s by searching for delta-functions in the second derivative of s
Parameters
----------
s : array_like
The signal
thresh : float
Threshold for delta search
delta : float
Sampling period
Returns
-------
intervals_s, peaks, dd : *intervals* are the intervals, *peaks* the found peaks and *dd* the second derivative of s.
"""
# derivative of derivative shows delta functions at discontinuities
dd = np.diff(np.diff(s, prepend=0), prepend=0)/delta**2
# find delta functions
peaks, _ = find_peaks(np.abs(dd), prominence=thresh)
# return interval limits in seconds
intervals = np.round(tau[peaks]*10)/10
return intervals, peaks, dd
def find_ind_least_diff(x, x0):
"""
Find index of the value that's nearest to x0
Parameters
----------
x : array_like
Array to be inspected
x0 : digit
x0 : float
Target value to be searched for
Returns
......@@ -32,7 +64,8 @@ def findIndOfLeastDiff(x, x0):
index : int or list
index(indices) of value(s) that is (are) nearest to x0
"""
if isinstance(x0, list) or isinstance(x0, np.ndarray):
return list(map(lambda t0: np.abs(x - t0).argmin(), x0))
else:
return np.abs(x - x0).argmin()
return np.abs(x - x0).argmin()
\ No newline at end of file
......@@ -153,7 +153,7 @@ def ient_ilaplace_ht(t=np.linspace(-6, 6, num=1024), H0=1, pp=np.array([]), pz=n
sd = np.array(a[0])
else:
# Zaehlergrad größer als Nennergrad. Keine Partialbruchzerlegung möglich
h = np.array([np.NaN for x in t])
h = np.ones(t.shape) * np.nan
return h, td, sd
......@@ -182,17 +182,20 @@ def ient_ilaplace_Hf(f=np.linspace(-6, 6, num=1024), H0=1, pp=np.array([]), pz=n
def ient_iz_ht(n=np.linspace(-6, 6, num=13), H0=1, pp=np.array([]), pz=np.array([]), ord_p=np.array([]),
ord_z=np.array([]), roc=[-12, 12]):
ord_z=np.array([]), roc=[0, 12]):
"""
Calculate inverse z-Transform from z-domain to n-domain
TODO - calculate according ord_p and ord_z, currently not used
"""
if np.isinf(roc[0]):
raise ValueError('Region of convergence must start at 0 or pole (=[0, b] anti-causal or =[b, inf] causal)')
# setting up parameters
poles = np.array(pp)
zeroes = np.array(pz)
#poles_order = np.array(ord_p)
#zeroes_order = np.array(ord_z)
nominator = H0 * np.poly(zeroes)
denominator = np.poly(poles)
......@@ -211,8 +214,9 @@ def ient_iz_ht(n=np.linspace(-6, 6, num=13), H0=1, pp=np.array([]), pz=np.array(
input_denominator = denominator
r, p, k = residuez(input_nominator, input_denominator)
#p[::-1].sort()
causal = p.real <= roc[0]
p = np.around(p, 5)
causal = np.around(np.abs(p), 5) <= np.abs(roc[0])
row = 0
......@@ -255,7 +259,7 @@ def ient_iz_Hf(f=np.linspace(-6, 6, num=1024), H0=1, pp=np.array([]), pz=np.arra
if dB:
return 20 * np.log10(np.maximum(eps, np.abs(numerator / denominator)))
else:
return np.abs(numerator / denominator), numerator, denominator
return np.abs(numerator / denominator)
def ient_ideal_sample(x, y, T):
def x_mirrored_around_zero(maxval, inc=1):
......