Pi (π) in Torch Tensor - torch.pi

There is no pi (π) function in PyTorch with some version, so some people suggested to use math.pi or np.pi for conversion into torch tensor. For me, torch.pi works with a notebook without the GPU, but not with a computer with the GPU.

Here is the code to compare different implementations of the Pi function.

import torch

import numpy as np

import time

import math

#test pi

time_start = time.perf_counter()

pi_math = torch.tensor(math.pi)

time_elapsed_test = (time.perf_counter() - time_start)

print('math pi time =',time_elapsed_test)

time_start = time.perf_counter()

pi_np = torch.tensor(np.pi)

time_elapsed_test = (time.perf_counter() - time_start)

print('np pi time =',time_elapsed_test)

time_start = time.perf_counter()

pi_torch = torch.pi

time_elapsed_test = (time.perf_counter() - time_start)

print('torch pi time =',time_elapsed_test)

time_start = time.perf_counter()

PI = torch.acos(torch.Tensor([-1]))

time_elapsed_test = (time.perf_counter() - time_start)

print('PI time =',time_elapsed_test)

Results for computation time:

math pi time = 3.6502000057225814e-05

np pi time = 1.9006000002264045e-05

torch pi time = 8.009999419300584e-07

PI time = 0.00035780000007434865

Hence torch.pi is the fastest.

Without torch.pi, you may use torch.tensor(np.pi)

For an accurate pi, cast the type to float64:

torch.tensor(np.pi, dtype=torch.float64)

Example:

(Pdb) torch.tensor(np.pi)*100000

tensor(314159.2812)

(Pdb) torch.tensor(np.pi, dtype=torch.float64)*100000

tensor(314159.2654, dtype=torch.float64)

MATLAB:

pi = 3.141592653589793

Reference:

Np.pi equivalent in Pytorch

Is there the Pi greek (3.1415…) defined somewhere?

Python: Short-Time Fourier Transform (STFT) and its inverse transform (ISTFT) using librosa

Below shows an example using STFT to transform speech into the frequency domain and ISTFT to transform the spectral data back to waveforms in the time domain. The Librosa package is used.

import librosa

import soundfile as sf

x, sr = librosa.load('1.wav', sr=None)

n_fft = 256

hop_length = 128

win_length = 256

X = librosa.stft(x, n_fft=n_fft, hop_length=hop_length, win_length=win_length)

y = librosa.istft(X, n_fft=n_fft, hop_length=hop_length, win_length=win_length)

sf.write('1_istft.wav', y, 16000)

Results for x and y:

References:

Matlab: Short-Time Fourier Transform (STFT) and its inverse transform (ISTFT) (StudyEECC)

librosa.stft

librosa.istft

Matlab: Short-Time Fourier Transform (STFT) and its inverse transform (ISTFT)

For sound processing in realtime scenarios, it is not possible to wait for a complete sound file. In this case, short-time processing is important.

Here is an example using STFT to transform speech into the frequency domain and inverse transform the spectral data back to waveforms in the time domain:

[S,F,T] = stft(x, fs,'Window',win2,'OverlapLength',128,'FFTLength',256); 

[y,ti] = istft(S,fs,'Window',win2,'OverlapLength',128,'FFTLength',256);

References:

stft - Short-time Fourier transform (MathWorks)

istft - Inverse short-time Fourier transform (MathWorks)

Matlab: Discrete Cosine Transform (DCT) for speech processing

The DCT and inverse DCT may be used to convert a speech signal into the transform domain using real values and back to the speech waveform.

Example with Matlab:

%load speech

load mtlb

x = mtlb;

X = dct(x);

y = idct(X);

Results:

The values for x and y look the same, but MATLAB consider them as different.

>> isequal(x,y)

ans =

  logical

   0

For real-time processing, short-time DCT is required.

References:

Discrete Cosine Transform (Wikipedia)

dct (MathWorks)

Matlab: Vector normalization to a new peak or RMS value

The following example shows how to normalize a vector to a new peak or RMS (root-mean-square) value. The rms function from the digital signal processing toolbox is used. For details about the calculations for RMS, see a previous article.

RMS = √(∑x_i²/n)

Vector a:

>> a = [1 3 10 2 5 7 8 17 14]

a =

     1     3    10     2     5     7     8    17    14

Normalize to a new peak:

>> target_peak = 15

target_peak =

    15

>> a_newpeak = a*target_peak/max(a)

a_newpeak =

    0.8824    2.6471    8.8235    1.7647    4.4118    6.1765    7.0588   15.0000   12.3529

Normalize to a new RMS value:

>> target_rms = 10

target_rms =

    10

>> a_newrms = a*target_rms/rms(a)

a_newrms =

    1.1051    3.3152   11.0506    2.2101    5.5253    7.7354    8.8405   18.7861   15.4709

Check the RMS value:

>> rms(a_newrms)

ans =

    10

References:

Root mean square (Wikipedia)/平方平均數(維基百科)

rms: Root-mean-square level (MathWorks)

Matlab: Min, Max, Mean, Square, Mean-Square, RMS (StudyEECC blog)

Matlab: Min, Max, Mean, Square, Mean-Square, RMS

The following example shows how to get the minimum, maximum, average/mean and RMS (root-mean-square) values from a vector.

RMS = √(∑x_i²/n)

>> a = [1 3 10 2 5 7 8 17 14]

a =

     1     3    10     2     5     7     8    17    14

>> b = min(a)

b =

     1

>> c = max(a)

c =

    17

>> d = mean(a)

d =

    7.4444

>> e = a.^2

e =

     1     9   100     4    25    49    64   289   196

>> f = mean(e)

f =


   81.8889

>> g = sqrt(f)

g =


    9.0492

where

b = minimum
c = maximum
d = mean/average
e = square
f = mean-square
g = root-mean-square (RMS)

The RMS value can also be worked out in a straight-forward command:

>> a_rms = sqrt(sum(a.^2)/length(a))

a_rms =

    9.0492

Or use rms in the digital signal processing toolbox:

>> a_rms_toolbox = rms(a)

a_rms_toolbox =

    9.0492

Reference:

Root mean square (Wikipedia)/平方平均數(維基百科)

rms: Root-mean-square level (MathWorks)

Matlab: Vector normalization to a new peak or RMS value (StudyEECC blog)

Matlab: Spectrogram of white noise processed by a 500-Hz Butterworth lowpass filter

Matlab code:

clc;clear;close all;

 

Fs = 16000;  %Sampling frequency

Fco = 500;   %Cutoff frequency

order = 4;

[b, a] = butter(order, Fco/(Fs/2), 'low');

 

%Generate white noise

x=2*rand(16000,1)-1;

 

y = filter(b, a, x); %get time domain result

 

spectrogram(x,[],[],[],Fs,'yaxis');colormap hot;title('x');figure

spectrogram(y,[],[],[],Fs,'yaxis');colormap hot;title('y');

Result:

References:

Divide a speech utterance into 3 bands

Apply a LPF in the time and frequency using filter() and freqz()

Matlab: Parseval's theorem

Parseval's theorem states that for N discrete points of signal,

Total Energy ∑ |x[n]|² = 1/N × ∑ |X[k]|²

where X[k] is the k^th point of the discrete Fourier transform of x[n].

Let's check this with Matlab:

load mtlb; %Load "Matlab" sound example

N = 128;

y = mtlb(1:N); %Read the 1st 128 points only

%Check Parseval's theorem

%y^2

y_2 = y.^2;

%Energy

energy = sum(y_2)

%Frequency domain

y_fft = fft(y); %Complex double

%|y_fft|^2

y_fft_2 = abs(y_fft).^2;

y_fft_2_check = y_fft.*conj(y_fft); %Should be identical to y_fft_2

energy_fft = sum(y_fft_2)/N

energy_fft_check = sum(y_fft_2_check)/N

Result:

Reference:

energy =

    2.5595


energy_fft =

    2.5595


energy_fft_check =

    2.5595

>> 

Parseval's theorem (Wikipedia)

Correlation Matrix 相關矩陣

correlation matrix 相關矩陣
一個穩態離散時間隨機過程可以一個時間序列來表示(M×1)，而當我們把這個隨機時間序列乘以它的赫密特轉置矩陣(1×M)，所得到的期望值(因為時間序列是隨機的)即是這個時間序列的相關矩陣(M×M)。

R = E[ u(n) u^H(n) ]

R: correlation matrix
u(n) : stochastic process
u^H(n) : Hermitian transpose
E: expectation

Hermitian transpose 赫密特轉置矩陣
將一個矩陣轉置並取共軛複數(complex conjugate)

在Matlab中，可直接以'符號求得Hermitian transpose，或使用ctranspose function

>> a = [1 1+j;2-j -2]

a =

   1.0000 + 0.0000i   1.0000 + 1.0000i
   2.0000 - 1.0000i  -2.0000 + 0.0000i

>> a'

ans =

   1.0000 + 0.0000i   2.0000 + 1.0000i
   1.0000 - 1.0000i  -2.0000 + 0.0000i


Hermitian matrix/self-adjoint matrix 埃爾米特矩陣/厄米特矩陣/自伴隨矩陣

共軛對稱的方塊矩陣，即 A = (A^T)^*，則A的Hermitian transpose等於A自己，即 A^H = A

Example:

>> a = [1 3+j; 3-j -2]

a =

   1.0000 + 0.0000i   3.0000 + 1.0000i
   3.0000 - 1.0000i  -2.0000 + 0.0000i

>> a'

ans =

   1.0000 + 0.0000i   3.0000 + 1.0000i
   3.0000 - 1.0000i  -2.0000 + 0.0000i

>> isequal(a,a')

ans =

  logical

   1

>> 

References:

Haykin, S. S. (2014). Adaptive filter theory. 5th edition, Pearson Education. pp 52-56.
Conjugate transpose (Wikipedia)
Hermitian matrix (Wikipedia) 埃爾米特矩陣(維基百科)

Adaptive Equalizer 自適應性等化器

在通訊系統中，例如對於電話數據機(modem)的資料傳輸，在低位元率(low bit rate)時，通常比較沒有符號間干擾/符際干擾(Intersymbol Interference, ISI)的問題。

但是當傳輸位元率提高時，例如超過2400 bits時，ISI的問題會變嚴重，這時就需要使用自適應性頻道等化器(Adaptive Channel Equalizer)來克服通道失真(Channel Distortion)的問題。

Adaptive Equalizer當見的實現方式是使用自適應性有限脈衝響應濾波器(Adaptive FIR Filter)，其係數採用LMS演算法(Least-Mean-Square Algorithm, 最小均方演算法)進行調整，以克服通道失真。LMS演算法是一種隨機梯度演算法(Stochastic Gradient Algorithm)。

通了使用LMS，Adaptive Equalizer也可採用RLS演算法(Recursive Least-Squares Algorithm, 遞迴最小平方演算法)。

參考資料：

Adaptive equalizer (Wikipedia)
Stochastic gradient descent (Wikipedia) (隨機梯度下降)
Least mean squares filter (Wikipedia) (屬於一種Stochastic Gradient Algorithm)
Recursive least squares filter (Wikipedia)

Haykin, S. S. (2014). Adaptive filter theory. 5th edition, Pearson Education. pp 40-43, 468-471.

Ingle, V. K., & Proakis, G. J. (2014). Essentials of Digital Signal Processing using MATLAB®. CENGAGE Learning, pp 603-606.

Matlab: FIR comb filter with various zero placement

The following Matlab code tests the effects of different number of zeros and checks the linear phase result.

a = 1;

b1 = [1 -1];

zplane(b1,a);
fvtool(b1,a);
LinPhase1 = islinphase(b1,a) %Check if linear phase

%Notch comb filter
b2 = [1 0 -1]; 
zplane(b2,a);
fvtool(b2,a);
LinPhase2 = islinphase(b2,a) %Check if linear phase

%Notch comb filter
b3 = [1 0 0 -1]; 
zplane(b3,a);
fvtool(b3,a);
LinPhase3 = islinphase(b3,a) %Check if linear phase

%Notch comb filter
b4 = [1 0 0 0 -1];
zplane(b4,a);
fvtool(b4,a);

%Notch comb filter
b5 = [1 0 0 0 0 -1];
zplane(b5,a);
fvtool(b5,a);
LinPhase5 = islinphase(b5,a) %Check if linear phase

Result:

LinPhase1 =

  logical
   1

LinPhase2 =

  logical
   1

LinPhase3 =

  logical
   1

LinPhase5 =

  logical
   1

b1 = [1 -1];

b2 = [1 0 -1]; 

b3 = [1 0 0 -1]; 

b4 = [1 0 0 0 -1]; 

 b5 = [1 0 0 0 -1]; 

Matlab: IIR Comb Notch or Peak Filter

IIR Comb Notch Filter

fs = 600;
fo = 100;
q = 35;
bw = (fo/(fs/2))/q;
[b,a] = iircomb(fs/fo,bw,'notch'); % Note type flag 'notch'
fvtool(b,a);
zplane(b,a);

Magnitude Response

Phase Response

Pole-Zero Diagram

IIR Comb Peak Filter

Replace the 'notch' string by 'peak':

[b,a] = iircomb(fs/fo,bw,'peak'); % Note type flag 'peak'

Magnitude Response

Phase Response

Pole-Zero Diagram

Reference

iircomb (MathWorks)

Matlab: Frequency response and pole-zero plot of a 1000-Hz Butterworth lowpass filter

The following Matlab code shows the magnitude and phase response and the Z-plane of a 1000-Hz Butterworth lowpass filter.

clc;clear;close all;

Fs = 8000;  %Sampling frequency
Fco = 1000; %Cutoff frequency
order = 4;
[b, a] = butter(order, Fco/(Fs/2), 'low');

zplane(b,a);title(['LPF Fco=',num2str(Fco),'Hz, Fs=',num2str(Fs),'Hz']);

w = 0:0.01:2*pi;

H = freqz(b, a, w);

mag = abs(H);
pha = angle(H);

figure
subplot(211);
plot(w,mag);title('Magnitude');xlabel('rad/s');
subplot(212)
plot(w,pha);title('Phase');xlabel('rad/s');

figure
w_Hz = w*Fs/(2*pi); %frequency in unit of Hz.
subplot(211);
plot(w_Hz,mag);title('Magnitude');xlabel('Hz');
subplot(212)
plot(w_Hz,pha);title('Phase');xlabel('Hz');

Result:

Pole-zero plot

Magnitude and phase responses with frequency in rad/s.

Magnitude and phase responses with frequency in Hz.

References:

How to plot the magnitude and phase of a given transfer function(z-domain)?
Pole-Zero plots of various filters
draw the Pole-Zero plot of the Z-transform

Illustration of Fundamental Frequency, Harmonics, and Formats

How to show the differences between the fundamental frequencies F0, harmonics and formants F1, F2, F3 etc. of a speech signal on a spectral diagram?

The following presentation file gives a good illustration.

http://research.cs.tamu.edu/prism/lectures/sp/l7.pdf

See slide #3.

Matlab: Pole-Zero plots of various filters

This page shows the comparisons of various Butterworth filter designs.

(1) Chang the sampling rate Fs = 4k, 8k, 16k. Fco = 200.

Fs = 2000 %4000, 8000, 16000
Fco = 200; %Cutoff frequency
order = 4;
[b, a] = butter(order, Fco/(Fs/2), 'low');

zplane(b,a);title(['LPF Fco=',num2str(Fco),'Hz, Fs=',num2str(Fs),'Hz']);

When the sampling frequency increases, the poles move towards to z = 1.

(2) Chang the cutoff frequency Fco = 200, 1k, 2k, 4k, 6k.

Fs = 16000
Fco = 200; %Cutoff frequency %200, 1000, 2000, 4000, 6000
order = 4;
[b, a] = butter(order, Fco/(Fs/2), 'low');

zplane(b,a);title(['LPF Fco=',num2str(Fco),'Hz, Fs=',num2str(Fs),'Hz']);

When the cutoff frequency increases, the poles move towards z = -1.

(3) Chang the order of the Butterworth filter. Order = 1, 2, 3, 4.

Fs = 16000
Fco = 1000;
order = 1; %1, 2, 3, 4
[b, a] = butter(order, Fco/(Fs/2), 'low');

zplane(b,a);title(['LPF Fco=',num2str(Fco),'Hz, Fs=',num2str(Fs),'Hz, Order = ',num2str(order)]);

The number of zeros = the order of the Butterworth filter.

(4) Chang the filter type: LPF, HPF, BPF, BSF.

Fs = 16000
Fco = 1000;
order = 4;
[b, a] = butter(order, Fco/(Fs/2), 'low'); %low, high

zplane(b,a);title(['LPF Fco=',num2str(Fco),'Hz, Fs=',num2str(Fs),'Hz']);

For the low pass filter, repeated zeros are at z=-1.

 For the high pass filter, repeated zeros are at z=1.

The band pass filter has repeated zeros at both z=1 and z=-1.

Fs = 16000
Fco1 = 1000; %3000
Fco2 = 2000; %4000
order = 4;
[b, a] = butter(order, [Fco1/(Fs/2),Fco2/(Fs/2)]);

zplane(b,a);title(['BPF Fco1=',num2str(Fco1),'Hz, Fco2=',num2str(Fco2),'Hz']);

The band stop filter (BSF) has repeated zeros on the unit circle. The zeros are surrounded by multiple poles.

Fs = 16000
Fco1 = 1000; %3000
Fco2 = 2000; %4000
order = 4;
[b, a] = butter(order, [Fco1/(Fs/2),Fco2/(Fs/2)],'stop');

zplane(b,a);title(['BSF Fco1=',num2str(Fco1),'Hz, Fco2=',num2str(Fco2),'Hz']);

Matlab: draw the Pole-Zero plot of the Z-transform

To plot the Pole-Zero diagram of the Z-transform on the Z-plane, write the transfer function in the form of b(z)/a(z) in terms of z^-n:

z(z-1)/(z+1/2)(z-1/4) = (1 - z^-1)/(1+1/2z^-1)(1-1/4 z^-1) =  (1 - z^-1)/(1+1/4 z^-1-1/8 z^-2)

Call the zplane() function to draw the ploe-zero plot with coefficients of b(z)/a(z):

>> b=[1,-1];
>> a=[1,1/4,-1/8];
>> zplane(b,a);

Result:

The zeros are at 0 and 1 and the poles are at -1/2 and 1/4.

Reference:

zplane (MathWorks)