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()

Compare the if condition between Matlab and Python

When coding with multiple programming languages, sometimes it is confusing to choose the right syntax for each language. Here is an example for the if condition in Matlab and Python:

Matlab

a = 1;

b = 2;

if a == 1

   disp('Hello');

else

   disp('Hi');

end

if b ~= 1

   disp('World'); 

end

Result:

Hello
World
>> 

Python

a = 1
b = 2

if a == 1:
   print('Hello')
else:
   print('Hi')

if b != 1:
   print('World')

Result:

Hello
World

STOI Code Test for Matlab and Python

(這篇文章主要是在測試Matlab和Python兩個版本程式所計算出的STOI值是否有差異)
STOI (Short-time objective intelligibility) is a prediction metric useful to replace previous objective prediction measures.

I tried to check if the Matlab and Python versions of STOI codes work equally with the identical output scores.

The Matlab code is available from Cees Taal's web site.

The python version is available from Pariente Manuel's GitHub.

The same processed wav file and reference wav file are tested using Matlab and Python.

The Python version needs a few configurations.

Create minicoda environments

conda create -n py38_STOI python=3.8 pytorch

conda activate py38_STOI

Install sound packages

pip install SoundFile

conda install -c conda-forge librosa

Run the python file.

python STOI_test.py

where the STOI_test.py is modified from Pariente Manuel's GitHub:

import soundfile as sf
from pystoi import stoi #or use from pystoi.stoi import stoi

clean, fs = sf.read('my_clean.wav')
denoised, fs = sf.read('my_processed.wav')

# Clean and den should have the same length, and be 1D
d = stoi(clean, denoised, fs, extended=False)
print(d)

And the result is

0.8157936197843473

The result of Matlab version is

>> d = stoi(y_my_clean, y_my_processed, fs)
d =
    0.8158

The results are similar, but there is still a tiny difference when showing more decimal points using digits() and vpa() commands.

>> digits(16)
>> vpa(d)
ans =
0.8157974711433416

The difference between Matlab and Python versions is small than 0.00001 and hence can be ignored.

References:

Cees Taal's web site
Pariente Manuel's GitHub

Matlab: Show more decimal points for a variable

To show more digits with Matlab, use digits() and vpa()

>> a = 1/3

a =

    0.3333

>> digits(6)
>> vpa(a)

ans =

0.333333

>> digits(10)
>> vpa(a)

ans =

0.3333333333

Reference:

digits (MathWorks)

Matlab: Log function with base 10

For the logarithmic function with base 10, simply use log10:

>> log10(100)

ans =

     2

>> 10*log10(0.1)

ans =

   -10

>> 20*log10(1000)

ans =

    60

>>

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)

Matlab: Paste figure into Microsoft Word or PowerPoint files

If Matlab figures are pasted into Microsoft Office applications such as Word and PowerPoint, the file becomes slow and the mouse icon is turning/loading.

I have figured out the following solution:

1. Select copy figure

2. Paste into an empty Word file.

3. Select the figure just pasted in the Word file and copy it.

4. Past into Gimp. You should NOT see this warning:

(If you see this error message, you probably directly pasted Matlab figure into Gimp. You cannot directly do this, so steps 2 and 3 are essential. The transparency can be preserved in this way.)

5. Select all in the image layer of Gimp and copy it.

6. Paste into the target Word/PowerPoint file.

Matlab: Singular and Nonsingular Matrices 奇異矩陣/非奇異矩陣

inverse matrix 逆矩陣/反矩陣

Nonsingular

If determinant ≠ 0 => i.e. the inverse matrix exists.
=> the matrix is invertible/non-singular.

invertible matrix 可逆矩陣
nonsingular matrix 非奇異矩陣/非特異矩陣

nondegenerate matrix 非退化矩陣

AB = BA = I

B = A^-1 = adj(A)/det(A)

where
adj(A) = adjoint matrix 伴隨矩陣 of A
det(A) = determinant 行列式 of A

>> a = [1 1 ; 2 1]

a =

     1     1
     2     1

>> det(a)

ans =

    -1

>> b = inv(a)

b =

    -1     1
     2    -1

>> a*b

ans =

     1     0
     0     1

>> b*a

ans =

     1     0
     0     1

Singular

If determinant = 0 => the matrix is invertible.

non-invertible matrix 不可逆矩陣

singular matrix 奇異矩陣/特異矩陣
degenerate matrix 退化矩陣

>> a = [2 3;1 1.5]

a =

    2.0000    3.0000
    1.0000    1.5000

>> det(a)

ans =

     0

>> b = inv(a)
Warning: Matrix is singular to working precision. 

b =

   Inf   Inf
   Inf   Inf

>>

Example:

References

奇異矩陣 singular matrix (國家教育研究院雙語詞彙、學術名詞暨辭書資訊網)
Invertible matrix (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) 埃爾米特矩陣(維基百科)

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: Find the index/row and column of the maximum element of a matrix

To find out the index of the maximum element of a matrix in Matlab, use:

>> myMat = [1 2 3 4;2 3 4 5;5 6 2 1]

myMat =
     1     2     3     4
     2     3     4     5
     5     6     2     1

>> [row, col] = find(myMat == max(myMat(:)))

row =
     3


col =
     2

>> myMat(3,2)

ans =
     6

>> max(myMat(:))

ans =
     6

Reference

How do I find the indices of the maximum (or minimum) value of my matrix? (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

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']);