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)
Information about Electrical, Electronic, Communication and Computer Engineering 電機、電子、通訊、電腦資訊工程的學習筆記
相關資訊~生醫工程:StudyBME
聽力科技相關資訊:電子耳資訊小站
iOS程式語言:Study Swift
樹莓派和Python:Study Raspberry Pi
2018/12/20
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
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.
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']);
(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)
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)
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)
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)
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) 邏輯迴歸 (維基百科)
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)
>> 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
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)
>> 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/07/04
Meta-Analysis 後設分析/整合分析
Meta-analysis (Wikipedia)
Meta-analysis (ScienceDirect)
整合分析 Meta-Analysis (國家教育研究院)
臨床醫師如何閱讀統合分析(Meta-analysis)的論文
童年經驗(1) 長大就好!? 42:20-43:15 (GoodTV 幸福學堂 錢玉芬)
Meta-analysis (ScienceDirect)
整合分析 Meta-Analysis (國家教育研究院)
臨床醫師如何閱讀統合分析(Meta-analysis)的論文
童年經驗(1) 長大就好!? 42:20-43:15 (GoodTV 幸福學堂 錢玉芬)
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)
References
Vector Space (Wikipedia)
What is a Vector Space? (Video)
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.
Result:
References:
stem - plot discrete sequence data (MathWorks)
Matlab: Unit Impulse δ[n] 單位脈衝函數
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)
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)
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)
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:
Result:
x =
1.0000 + 2.0000i
ans =
-2.0000 + 1.0000i
ans =
2.0000 - 1.0000i
Related Information:
Complex numbers - simple calculations (StudySwift - iOS)
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/維基百科)
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.
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)
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
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/13
Matlab: Differentiation of sine and cosine functions
To calculate the derivates of sin() and cos(), use the syms symbolic variable as below:
>> syms x
>> y = sin(x);
>> y_prime = diff(y)
y_prime =
cos(x)
>> y_prime_prime = diff(y_prime)
y_prime_prime =
-sin(x)
References
Matlab: Basic Differentiation and Integration with the syms symbolic variable
Basic Calculus Terms 基本微基分名詞
>> syms x
>> y = sin(x);
>> y_prime = diff(y)
y_prime =
cos(x)
>> y_prime_prime = diff(y_prime)
y_prime_prime =
-sin(x)
References
Matlab: Basic Differentiation and Integration with the syms symbolic variable
Basic Calculus Terms 基本微基分名詞
Basic Calculus Terms 基本微基分名詞
calculus 微積分
differentiation 微分
integration 積分
derivative 導數
integral 積分
References
Glossary of calculus (Wikipedia)
微積分常用述語中英文對照
Matlab: Basic Differentiation and Integration with the syms symbolic variable
Matlab: Differentiation of sine and cosine functions
differentiation 微分
integration 積分
derivative 導數
integral 積分
References
Glossary of calculus (Wikipedia)
微積分常用述語中英文對照
Matlab: Basic Differentiation and Integration with the syms symbolic variable
Matlab: Differentiation of sine and cosine functions
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.
References
Glossary of calculus (Wikipedia)
微積分常用述語中英文對照 (交大微積分教學小組)
Basic Calculus Terms 基本微基分名詞 (Study EECC)
Matlab: Differentiation of sine and cosine functions
>> 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.
References
Glossary of calculus (Wikipedia)
微積分常用述語中英文對照 (交大微積分教學小組)
Basic Calculus Terms 基本微基分名詞 (Study EECC)
Matlab: Differentiation of sine and cosine functions
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"
>> 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)
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)
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)
C++ (helloCPP.cpp)
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
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");
}
#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)
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/)
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
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
>> 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
訂閱:
文章 (Atom)