2018/12/20

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)

2018/12/15

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

2018/11/29

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.

2018/11/06

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)

2018/10/06

Matlab: Sinc function

The code below shows how to use the sinc() function of Matlab's Signal Processing Toolbox:

x = (-6:1/100:6);
y1 = sinc(x);
plot(x,y1);
y2 = sinc(2*x);
plot(x,y2);
xlabel('x');ylabel('sinc'); title('Sinc Function')
legend('sinx(x)', 'sinc(2x)')

Result:

Note:
The sinc() function here is defined as

    sinc(x) = sin(πx)/πx (normalized sinc function)

References:

sinc - sinc function (MathWorks)
MATLAB sinc issue (StackOverflow)

2018/09/04

Matlab: Find the arguments of maxima/minima (argmax/argmin) of a vector

The argument of maximum means the index at which the vector/function values is maximized.

Example

>> a = [3 2 1 -2 7]

a =
     3     2     1    -2     7

>> [max_value, argmax] = max(a)

max_value =
     7

argmax =
     5

>> [min_value, argmin] = min(a)

min_value =
    -2

argmin =
     4

>> 

References:

Arg max (Wikipedia)
how do i find argmax? (MathWorks)

2018/09/03

Matlab: feval() - Call functions with the function name as a string input

To call any Matlab function, use feval() with the function nameas the first input string and the parameters of the function as the other remaining input arguments.

For example, for sayHello.m:

function sayHello(name)

fprintf('Hello %s!\n',name);

Results of the command line with fevel():

>> sayHello('Matlab')
Hello Matlab!
>> feval('sayHello','World')
Hello World!

Reference:

feval() - Evaluation function (MathWorks)

2018/08/21

Mathematics for Machine Learning 機器學習的數學

Siraj Raval的影片介紹,機器學習主要有以下幾種的數學:

Calculus 微積分 - 最佳化

Linear Algebra 線性代數 - 實現演算法

Probability 機率 - 預估結果

Statistics 統計 - 找出目標

----

參考資料

Mathematics of Machine Learning (Siraj Raval)
Logistic regression (Wikipedia) 邏輯迴歸 (維基百科)

2018/08/09

Matlab: Find local maxima 求局部極大值

To find the local maximum values in Matlab, use findpeaks():

>> a = [7 -2 -5 8 -6 -3 4 9 6 3 -2 5 8 11 7 3 -5];
[peaks, items] = findpeaks(a);
>> peaks

peaks =

     8     9    11

>> items

items =

     4     8    14

>> 

To find N large elements of a vector in Matlab in descending order, use maxK():

K>> [peaks, items] = maxk(a,5);
K>> peaks

peaks =

    11     9     8     8     7

K>> items

items =

    14     8     4    13     1

K>> 

Reference:

findpeaks - Find local maxima (MathWorks)

2018/08/07

Matlab: Flip vector 水平翻轉向量

>> a = [1 2 3 4 5 6];
>> flip(a)

ans =

     6     5     4     3     2     1

2018/07/30

Matlab: Maximum element and index

To the the maximum value of an array and the index of the maximum value, use the max() function  as below:

>> a = [1 4 5 7 3 2]

a =

     1     4     5     7     3     2

>> [max, index] = max(a)

max =

     7


index =

     4

>> 

Reference:

max - Maximum elements of an array (MathWorks)

2018/06/28

Vector Spaces 向量空間

2D Vector Space
vector v = (x, y) or (v1, v2)

3D Vector Space
vector v = (x, y, z) or (v1, v2, v3)

N-Dimensional Vector Space
vector v = (v1, v2, v3, ..., vN)

References

Vector Space (Wikipedia)
What is a Vector Space? (Video)

2018/06/16

Matlab: Plot discrete sequence data with STEM() 繪製離散序列資料

To plot discrete sequence data in Matlab, use the stem() function:

Example:
For y = e-x + 1 with x between -2 and 2, plot the continuous y curve and discrete stem diagram.

x = -2:0.1:2;
y = exp(-x)+1;
subplot(211);
plot(x,y);
xlabel('x');ylabel('y');title('y = exp(-x)+1');
subplot(212);
stem(x,y);

xlabel('x');ylabel('y');title('discrete x');

Result:



References:

stem - plot discrete sequence data (MathWorks)
Matlab: Unit Impulse δ[n] 單位脈衝函數

2018/06/15

Matlab: Resample 重新取樣

To change the sample rate of an audio signal with Matlab, use the resample function:

load mtlb;
sound(mtlb,Fs);
pause(1); %Wait until the sound is finished
Fs_new = 16e3; %change the sample rate to 16000
[P,Q] = rat(Fs_new/Fs);

y = resample(mtlb,P,Q);

sound(y,Fs_new);

Reference

Changing Signal Sample Rate (MathWorks)

2018/06/11

Floating-Point

IEEE 754-1985, 754-2008
IEEE single-precision standard (binary32) - 32 bits
IEEE double-precision standard (binary64) - 64 bits

single-recision (binary32) -
1 bit - sign (S)
8 bits - exponention (E)
23 bits - mantissa (M)

Number = (-1)S×1.M×E-127

Reference:

Tutorial: Floating-Point Binary
Floating point ALU using VHDL implemented on FPGA
Single-precision floating-point format (Wikipedia)
Double-precision floating-point format (Wikipedia)

Two's complement

For representing +ve and -ve numbers.

The most significant bit (leftmost) bit:
0 => 0 or positive
1 => negative

e.g. 01110000 => +ve
e.g. 11110000 => -ve

If the left most bit = 0,

01110000 = 2^6+2^5+2^4 = 64+32+16 = 112

If the left most bit = 1, invert every remaining bit and add 1.

11110000 => 0001111 + 1 => 0010000 => -16

Reference:

Two's complement (Wikipedia)

2018/06/07

Matlab: Real and imaginary parts of complex number and multiplication with imaginary unit j

For a complex number x = 1 + 2i or x = 1 + 2j

Matlab commands for real part - real(x)
Matlab commands for imaginary part - imag(x)

Example:

>> x = 1 + 2j

x =

   1.0000 + 2.0000i

>> real(x)

ans =

     1

>> imag(x)

ans =

     2

>>

To multiply x with an imaginary unit i or -i, simply swap the real part with the imaginary part and add a negative sign.

Example:


x = 1 + 2j
mult_pos_j(x)
mult_neg_j(x)

function result = mult_pos_j(x)
    % x = a + bj
    % x(j) = aj - b = -b + aj
    result = -imag(x) + real(x)*j;
end

function result = mult_neg_j(x)
    % x = a + bj
    % x(-j) = -aj + b = b - aj
    result = imag(x) - real(x)*j;
end

Result:

x =

   1.0000 + 2.0000i


ans =

  -2.0000 + 1.0000i


ans =

   2.0000 - 1.0000i

Related Information:

Complex numbers - simple calculations (StudySwift - iOS)

2018/05/31

Identity Matrix 單位矩陣

Identity Matrix / Unit Matrix 單位矩陣/恆等矩陣

Example:

I =
[ 1 0 ]
[ 0 1 ]

I =
[ 1 0 0 ]
[ 0 1 0 ]
[ 0 0 1 ]

I =
[ 1 0 0 0 ]
[ 0 1 0 0 ]
[ 0 0 1 0 ]
[ 0 0 0 1 ]

Matlab Example:

>> I = eye(4)

I =

     1     0     0     0
     0     1     0     0
     0     0     1     0
     0     0     0     1

>>

For an NxN identity matrix I and a NxN matrix A,

AI = IA = A

For an NxN identity matrix I and a Nx1 matrix B or a 1xN matrix C,

IB = B
CI = C

Reference:

Identity Matrix (Wikipedia/維基百科)

2018/05/05

Matlab: cross-correlation and autocorrelation 互相關與自相關

cross-correlation 互相關
autocorrelation 自相關

For autocorrelation function in Matlab:
autocorr - Econometrics Toolbox - for statistics. 😳
xcorr - Signal Processing Toolbox - for engineers!! 😄

To use xcorr:
xcorr(x, y) => cross-correlation of discrete signals x and y.
xcorr(x) => autocorrelation of x. Special condition of cross-correlation since auto correlation is the cross-correlation of signal x and shifted version of itself.

Autocorrelation is a special case of cross-correlation.

Autocorrelation may be used for pitch (fundamental frequency, F0) detection.


References:

Cross-correlation (Wikipedia)
Autocorrelation (Wikipedia)
Pitch detection algorithm (Wikipedia)
D. Gerhard. Pitch Extraction and Fundamental Frequency: History and Current Techniques, technical report, Dept. of Computer Science, University of Regina, 2003.

2018/04/18

Signal-to-Noise Ratio (SNR)

Signal-to-Noise Ratio (SNR) 訊號雜訊比/訊雜比/訊噪比/信噪比

SNR = (signal power)/(noise power) = Psignal/Pnoise
SNRdB = 10log10(SNR) = 10log10(Psignal)-10log10(Pnoise) in decibels (dB)

Therefore, if Psignal > Pnoise,

SNR > 1 and SNRdB > 0 dB.

If SNRdB < 0, Psignal < Pnoise. In other words, when the SNRdB is negative, the power of the signal is weaker than the power of background noise.

Reference

Signal-to-noise ratio/訊號雜訊比 (Wikipedia)

2018/04/14

Matlab: Laplace Transform

The examples below demonstrate the Laplace transforms and inverse Laplace transforms using Matlab symbolic number or variable.

Laplace Transform

L{1} = 1/s

>> x = sym(1);
>> laplace(x)

ans =

1/s

L{t} = 1/s2

>> syms t
>> laplace(t)

ans =

1/s^2

Inverse Laplace Transform

L-1{1/s} = 1

>> syms s
>> ilaplace(1/s)

ans =

1

L-1{1/(s-a)} = eat where a = 2

>> ilaplace(1/(s-2))

ans =

exp(2*t)

L-1{a2/(s2+a2)} = sin(at) where a = 1

>> ilaplace(1/(s^2+1))

ans =

sin(t)

Reference

Zill, D., Wright, W. S., & Cullen, M. R. (2011). Advanced engineering mathematics. 4th Edition. Jones & Bartlett Learning. p 197-198.

laplace / ilaplace - MathWorks Documentation

2018/04/12

Matlab: Basic Differentiation and Integration with the syms symbolic variable

To perform basic calculus calculations with Matlab, declare symbolic variables with syms command:

>> syms x;
>> y = 3*x^3+2*x^2-4*x+5;

Hence y = 3x3+2x2-4x+5.

For differentiation, the derivative of y is y' = diff(y):

>> y_prime = diff(y)

y_prime =

9*x^2 + 4*x - 4

Hence y' = dy/dx = 9x2+4x-4.

For integration, the integral of y is ∫ y dy = int(y):

>> y_integral = int(y)

y_integral =

(3*x^4)/4 + (2*x^3)/3 - 2*x^2 + 5*x

Consequently, ∫ y dy = (3x4)/4 + (2x3)/3 - 2x2 + 5x + c.

2018/04/10

Matlab: Simple structure with dot (.)

To create a basic form of structure in Matlab, simply use a dot (.) as below:

>> mystruct.name = "peter";
>> mystruct.age = 15;
>> mystruct.hobby = "swimming";
>> mystruct

mystruct =

  struct with fields:

     name: "peter"
      age: 15
    hobby: "swimming"




2018/03/29

Verilog Development Environment with Mac

To develop Verilog codes with a Mac, the following tools may be used:

Command line editor: Nano

Graphical UI Editor: TextWrangler (Available on App Store for Free)

Compiler: Icarus Verilog (Install with Homebrew with command: brew install icarus-verilog)

Waveform Viewer: GTKWave (Write $dumpfile command to output vcd file

References:

如何在Mac OS X上安裝Verilog環境 (Instructions to Icarus Verilog and GTKWave installation)
is there a good verilog editor for mac?

Verilog vs. VHDL (Study EECC)
Digital Logic Design - Truth Table and K-map (Study EECC)

2018/03/16

Digital Logic Design - Truth Table and K-map

To convert digital logic conditions in a truth table to a digital circuit, use the K-map.

truth table 真值表
Karnaugh map / K-map 卡諾圖
Boolean algebra / logic algebra 布林代數 / 邏輯代數

References

Karnaugh map (Wikipedia)
Truth table (Wikipedia)
Boolean algebra (Wikipedia)
HOW TO: Combinational logic: Truth Table → Karnaugh Map → Minimal Form → Gate Diagram (YouTube)

Verilog Development Environment with Mac (Study EECC)

2018/03/10

Matlab: Run Hello World files written in C/C++ with MEX command

To run C/C++ files from Matlab, use the MEX command.
Create the files in C or C++:

C (hello.c)

#include"mex.h"
#include<stdio.h>

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
    printf("Hello, World!\n");
}


C++ (helloCPP.cpp)

#include"mex.h"
#include <iostream>

using namespace std;

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
    cout << "Hello, World!\n";

}

The mexFunction in either of the C or C++ example above is similar to the main function in pure C or C++ programs. Type your main program here such as printf in C or cout in C++.

Execute the files with the MEX command in Matlab such as below:

C
mex hello.c
hello

C++
mex helloCPP.cpp
helloCPP

Results in the Matlab command window:


References

Matlab 教材:測試 MEX-file
C language with Mac: Hello World with Mac's Terminal
C++ with Mac: Hello World with Mac's Terminal

2018/02/23

Verilog vs. VHDL

Verilog -
1984年推出、1995年成為標準IEEE 1364-1995、較廣泛使用、接近程式語言

VHDL (VHSIC Hardware Description Language) -
1983開始、1987年成為標準IEEE 1076-1995、使用較Verilog少、學習時間較Verilog長、仍有其優點

參考資料
Verilog (Wikipedia)
VHDL (Wikipedia)
SystemVerilog (Wikipedia)
Verilog HDL和VHDL的比較

Verilog Development Environment with Mac (Study EECC)

2018/02/06

Matlab: Generate a pure tone with *.m file

To play a pure tone (sine wave) in Matlab, type the following code in a new file and save it as genTone.m:

function genTone(Fs, F, time)
t = 0:1/Fs:time;
y = sin(2*pi*F*t);
sound(y,Fs);
end

Now play sound in command

>> genTone(16000,100,0.5);
>> genTone(16000,1000,0.5);

The 3000Hz sound below will be wrong since it is greater than the 4000/2 = 2000 Hz Nyquist frequency.

>> genTone(4000,3000,0.5);

References:

how do i generate sound using MATLAB?
Chapter 2: Practical Audio Processing / tonegen.m (http://mcloughlin.eu/)

2018/01/26

Matlab: Divide a speech utterance into 3 bands


>> load mtlb %Load the 'Matlab' speech / dimention: 4001x1
>> sound(mtlb,Fs); %Play the original sound
>> len = length(mtlb); %Number of points
>> time = [1:len]/Fs; %x-axis time vector / dimention: 1x4001
>> subplot(5,1,1);
>> plot(time,mtlb);
>> title('Original Speech');

>> Fco1 = 1000/(Fs/2); %Cutoff freq. 1
>> Fco2 = 2000/(Fs/2); %Cutoff freq. 2

>> [Bhpf, Ahpf] = butter(4,Fco2,'high'); %design a HPF
>> [Bbpf, Abpf] = butter(4,[Fco1 Fco2],'bandpass'); %design a BPF
>> [Blpf, Alpf] = butter(4,Fco1,'low'); %design a LPF


>> y_hpf = filter(Bhpf, Ahpf, mtlb); %get time domain result

>> y_bpf = filter(Bbpf, Abpf, mtlb); %get time domain result
>> y_lpf = filter(Blpf, Alpf, mtlb); %get time domain result

>> subplot(5,1,2);
>> y_reconstruct = y_hpf + y_bpf + y_lpf;
>> plot(time,y_reconstruct);
>> title('Reconstructed Speech');
>> sound(y_reconstruct,Fs); %Play the reconstructed sound

>> subplot(5,1,3);
>> plot(time,y_hpf);
>> title('HPF Speech');
>> sound(y_hpf,Fs); %Play the HPF sound


>> subplot(5,1,4);
>> plot(time,y_bpf);
>> title('BPF Speech');
>> sound(y_bpf,Fs); %Play the BPF sound


>> subplot(5,1,5);
>> plot(time,y_lpf);
>> title('LPF Speech');
>> sound(y_lpf,Fs); %Play the LPF sound


Result:


2018/01/06

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

The code below shows how to implement a 1000Hz Butterworth low pass filter using the filter() function in time domain and freqz() function in the complex frequency domain.

>> load mtlb %Load the 'Matlab' speech / dimention: 4001x1
>> sound(mtlb,Fs); %Play the original sound
>> len = length(mtlb); %Number of points
>> time = [1:len]/Fs; %x-axis time vector / dimention: 1x4001
>> subplot(3,1,1);
>> plot(time,mtlb);
>> title('Original Speech');

>> [b, a] = butter(3,1000/(Fs/2),'low'); %design a Butterworth filter


Time domain

>> y1 = filter(b, a, mtlb); %get time domain result directly % dimention: 4001x1
>> sound(y1,Fs); %Play the sound processed by filter()
>> subplot(3,1,2);
>> plot(time,y1);
>> title('filter()');

Frequency domain (Method 1)
>> mtlb_FFT = fft(mtlb); %transformed input signal / dimention: 4001x1
>> H = freqz(b, a, len, 'whole'); %filter transfer function / dimention: 4001x1
>> Y2_FFT = mtlb_FFT .* H; %filtered result in freq. domain = multiplication of transformed speech signal and the filter transfer function in the freq. domain. % dimention: 4001x1
>> y2 = ifft(Y2_FFT); %Inverse FFT to time domain % dimention: 4001x1
>> sound(y2, Fs);
>> subplot(3,1,3);
>> plot(time,y2);
>> title('freqz() and transform');

Frequency domain (Method 2)
>> mtlb_FFT = fft(mtlb); %transformed input signal / dimention: 4001x1
>> H = freqz(b, a, len); %filter transfer function / dimention: 4001x1
>> Y2_FFT = mtlb_FFT .* H; %filtered result in freq. domain = multiplication of transformed speech signal and the filter transfer function in the freq. domain. % dimention: 4001x1
>> y2 = ifft(Y2_FFT); %Inverse FFT to time domain % dimention: 4001x1
>> sound(abs(y2), Fs);
>> subplot(3,1,3);
>> plot(time,y2);
>> title('freqz() and transform');

Result:


References:

Matlab: Low pass filtered Chirp sound - Study EECC
Filter Applications (濾波器應用) - Audio Signal Processing and Recognition (音訊處理與辨識)
MATLAB: filter in the frequency domain using FFT/IFFT with an IIR filter - StackOverflow