🔀 switch the default branch to 'main' from 'master'

2016 Fall/CommE5030 - Time Frequency Analysis and Wavelet Transform/README.md

@@ -0,0 +1,54 @@
+# Course Information
+
+<p align="center">
+  <img src="https://i.imgur.com/EM0m4TD.png">
+</p>
+
+## Course Name
+
+時頻分析和小波轉換
+Time-frequency Analysis and Wavelet Transform
+
+## Semester
+
+105-1
+
+## Course Number
+
+CommE 5030
+
+## Description
+
+時頻分析可視為傅立葉轉換的一般化，本課程將介紹多個時頻分析的方法，包括短時傅立葉轉換，加伯轉換，韋格納分佈，柯恩轉換，希爾伯特黃轉換，小波轉換，並介紹它們在濾波器設計，聲學，影像處理，通訊，和壓縮感知上的應用。
+
+## Objective
+
+一方面介紹時頻分析的理論和應用，一方面也複習信號處理的知識。
+
+## Requirement
+
+- **平時分數**：15 Scores  
+    根據上課回答問題加分
+- **作業分數**：60 Scores  
+    共 5 次, 每 3 週一次
+- **報告分數**：25 Scores  
+    方式有四種
+    1. 書面報告（10 頁以上，不含封面），中英文皆可，使用 11 或 12 的字體，題目可選擇和課程有關的任何一個主題，格式和一般寫期刊論文或碩博士論文相同，包括 Abstract, Conclusion, 及 References，並且要分 Sections，必要時有 Subsections。
+    2. Tutorial（和書面報告格式相同，但 18 頁以上，題目由老師指定，以清楚的介紹一個主題的基本概念和應用為要求，選擇這個項目的同學，學期成績加 3 分）
+    3. 口頭報告（限四個人，每個人 30~40 分鐘，題目可選擇和課程有關的任何一個主題，選擇這個項目的同學，學期成績加 2 分）
+    4. 編輯 Wikipedia（中文或英文網頁皆可，至少 2 個條目，但不可同一條目翻成中文和英文。限和課程相關者，自由發揮，越有條理、有系統的越好）
+
+## Reference
+
+### 授課用書：
+
+- 上課講義
+
+### 參考書籍： 
+
+1. S. Qian, D. Chen, Joint Time-Frequency Analysis: Methods and Applications, Prentice-Hall, 1996.
+2. S. Mallat, A Wavelet Tour of Signal Processing: The Sparse Way, Academic Press, 3rd edition, 2009.
+3. L. Cohen, Time-Frequency Analysis, Prentice-Hall, New York, 1995.
+4. K. Grochenig, Foundations of Time-Frequency Analysis, Birkhauser, Boston, 2001.
+5. L. Debnath, Wavelet Transforms and Time-Frequency Signal Analysis, Birkhäuser, Boston, 2001.
+

2016 Fall/CommE5030 - Time Frequency Analysis and Wavelet Transform/homework/homework_01/hw01.m

@@ -0,0 +1,68 @@
+function hw01()
+defsel = input('Do you want to define the variable and input x by yourself? (Y/N)','s');
+switch defsel
+    case 'Y'
+        dt = input('Input the Delta t: ');
+        df = input('Input the Delta f: ');
+        freq = input('Input the frequency interval: ');
+        t1 = input('Input the t2 (Time for signal 1): ');
+        t2 = input('Input the t2 (Time for signal 2): ');
+        t3 = input('Input the t3 (Time for signal 3): ');
+        t = input('What is the samples on t-axis?: ');
+        f = input('What is the samples on f-axis?: ');
+        x = input('What is the input: ');
+        B = input('What is the interval [-B,B] of integration?: ');
+    case 'N'
+        dt = 0.05;
+        df = 0.05;
+        freq = [-1500:1500];
+        t1 = [0:dt:10-dt];
+        t2 = [10:dt:20-dt];
+        t3 = [20:dt:30];
+        t = [0:dt:30];
+        f = [-5:df:5];
+        x = [cos(2*pi*t1),cos(6*pi*t2),cos(4*pi*t3)];
+        B = 1;
+end
+tic
+y = recSTFT(x,t,f,B);
+toc
+drawgraph(y,freq,t)
+end
+
+function [y] = recSTFT ( x,t,f,B )      % Define recSTFT function
+T = length(t);                          % Get the length of vector T
+F = length(f);                          % Get the length of Vector F
+delta_f = f(2) - f(1);                  % Calculate Delta f
+delta_t = t(2) - t(1);                  % Calculate Delta T
+Q = B/delta_t;
+x = [ x, zeros(1,Q) ];                  % Add zero
+N = 1 / (delta_f * delta_t);
+X = zeros(T,F);
+m = f / delta_f;
+n = t / delta_t;
+m1 = mod(m,N) + 1;
+for i = 1:T
+    x1 = zeros(1,floor(N));
+    if i <= Q+1
+        x1((Q-i+2):(2*Q+1)) = x(1:(Q+i));
+    else
+        x1(1:(2*Q+1)) = x((i-Q):(i+Q));
+    end
+    X1 = fft(x1);
+    XX1 = zeros(1,length(m));
+    XX1 = X1(floor(m1));
+    X(i,:) = XX1*delta_t.*exp(1j*2*pi*(Q-(i-1))*m/N);
+end
+y = X';
+end
+
+function drawgraph(y,freq,t)
+image(t,freq,abs(y)/max(max(abs(y)))*400)      % Here 400 is a constant can be changed.
+colormap(gray(256))                     % Take color in winter set.
+set(gca,'Ydir','normal')                % Make the y-axis upsidedown.
+set(gca,'Fontsize',12)                  % Change the font size.
+xlabel('Time (Sec)','Fontsize',12)      % Define x-axis.
+ylabel('Frequency (Hz)','Fontsize',12)  % Define y-axis.
+title('STFT of x(t)','Fontsize',12)     % Define the graph title.
+end

2016 Fall/CommE5030 - Time Frequency Analysis and Wavelet Transform/homework/homework_02/hw02.m

@@ -0,0 +1,63 @@
+function hw02()
+funsel = input('Function wavread for Older Matlab / audioread for Matlab 2016. (W/A)','s');
+switch funsel
+    case 'W'
+		[a1, fs] = wavread('Chord.wav');		% read audio file
+    case 'A'
+		[a1, fs] = audioread('Chord.wav');		% read audio file
+end
+x = a1(:,1);
+tau = 0:1/44100:1.6077;
+dt = 0.01;
+df = 1;
+t = 0:dt:max(tau);
+f = 20:df:1000;
+sgm = 200;
+tic												% Timer Start
+y = Gabor(x, tau, t, f, sgm);					% Define gabor function
+toc												% Timer End
+drawgraph(y,f,t)
+end
+
+function y = Gabor( x, tau, t, f, sgm )
+T = length(tau);
+F = length(f);
+C = length(t);
+dt = t(2) - t(1);
+df = f(2) - f(1);
+dtau = tau(2) - tau(1);
+N = 1/(df * dtau);
+B = 1.9143/(sgm^(1/2));							% |a| > 1.9143, the Gaussian Function is small enough to be ignored.
+Q = round(B / dtau);
+n0 = tau(1) / dtau;
+m0 = f(1) / df;
+c0 = t(1) / dt;
+S = dt / dtau;
+y = zeros(C,F);
+x1 = zeros(1,N);
+window = (sgm^(1/4)).* exp(-sgm.*pi.*((Q-(0:N-1))*dtau).^2);
+for n = c0:(c0+C-1)
+    win_const = dtau*exp( (1j*2*pi*(Q-n*S)).*(m0:(m0+F-1))./N);
+    ns_q = n*S-Q;								% Define variable ns_q = (n*S)-Q
+    for q = 0:N-1
+        if (q <= (2*Q) && (ns_q+q>=0) && ((ns_q+q+1)<=T))
+            x1(q+1) = x(ns_q+q+1);
+        else
+            x1(q+1) = 0;
+        end
+    end
+    fft_out = fft((window.*x1),N);
+    y(n+1-c0,:) = ( win_const .* fft_out(m0+1:m0+F))';
+end
+y=y';
+end
+
+function drawgraph(y,f,t)
+image(t,f,abs(y)/max(max(abs(y)))*400)			% Here 400 is a constant can be changed.
+colormap(gray(256))								% Take color in gray.
+set(gca,'Ydir','normal')						% Make the y-axis upsidedown.
+set(gca,'Fontsize',12)							% Change the font size.
+xlabel('Time (Sec)','Fontsize',12)				% Define x-axis.
+ylabel('Frequency (Hz)','Fontsize',12)			% Define y-axis.
+title('Gabor Transform','Fontsize',12)			% Define the graph title.
+end

2016 Fall/CommE5030 - Time Frequency Analysis and Wavelet Transform/homework/homework_03/hw03.m

@@ -0,0 +1,45 @@
+function hw03()
+% Define Start
+t = -9:0.0125:9;
+f = -4:0.05:4;
+x = exp(j*t.^2/10-j*3*t).*((t>=-9)&(t<=1))+exp(j*t.^2/2+j*6*t).*exp(-(t-4).^2/10);
+tic                 % Timer Start
+y = wdf(x, t, f);   % Define gabor function
+toc                 % Timer End
+drawgraph(y,f,t)
+end
+
+function y = wdf(x,t,f)
+dt = t(2)-t(1);
+df = f(2)-f(1);
+N = round(1/(2*dt*df));
+n1 = round(t(1)/dt);
+n2 = round(t(length(t))/dt);
+m1 = round(f(1)/df);
+m2 = round(f(length(f))/df);
+m=mod([m1:m2],N)+1;
+Lt = n2-n1+1;
+Lf = m2-m1+1;
+y = zeros(Lf,Lt);
+for n = n1:n2
+    U = min(n2-n,n-n1);
+    Q = 2*U+1;
+    A = x(1-n1+[n-U:n+U]).*(x(1-n1+[n+U:-1:n-U])').';
+    A1 = fft(A,N)*2*dt;
+    a1 = ceil(Q/N);
+    for a2 = 2:a1
+        A1 = A1+fft(A((a2-1)*N+1:min(a2*N,Q)),N)*2*dt;
+    end
+    y(:,n-n1+1)=(A1(m).*exp(j*2*pi/N*U*(m-1))).';
+end
+end
+
+function drawgraph(y,f,t)
+image(t,f,abs(y)/max(max(abs(y)))*400)						% Here 400 is a constant can be changed.
+colormap(gray(256))											% Take color in gray.
+set(gca,'Ydir','normal')									% Make the y-axis upsidedown.
+set(gca,'Fontsize',12)										% Change the font size.
+xlabel('Time (Sec)','Fontsize',12)							% Define x-axis.
+ylabel('Frequency (Hz)','Fontsize',12)						% Define y-axis.
+title('Wigner Distribution Function','Fontsize',12)			% Define the graph title.
+end

2016 Fall/CommE5030 - Time Frequency Analysis and Wavelet Transform/homework/homework_04/hw04.m

@@ -0,0 +1,63 @@
+function hw04()
+t = [0: 0.01: 10];
+x = 0.2*t + cos(2*pi*t) + 0.4*cos(10*pi*t);
+thr = 0.2;
+tic													% Timer Start
+y = hht(x, t, thr);
+toc													% Timer End
+drawgraph(y, t)
+end
+
+function y = hht(x, t, thr)
+n = 1;
+dt = t(2) - t(1);
+while (1)
+    if (length(findpeaks(x)) <= 3)
+        y(n,:) = x;
+        break
+    end
+
+    temp = x;
+    test = 1;
+    k = 1;
+    while (test == 1 && k < 3)
+        test = 0;
+        [max maxloc] = findpeaks(temp);         	% Step02
+        peaks = spline((maxloc-1)*dt, max, t);    	% Step03
+        [neg_min minloc] = findpeaks(temp*(-1));	% Step04
+        min = (-1)*(neg_min);
+        dips = spline((minloc-1)*dt, min, t);     	% Step05
+        z = (peaks + dips) / 2;                    	% Step06_1
+        h = (temp - z);								% Step06_2
+
+        % Step07
+        hpeaks = findpeaks(h);
+        hdips = (-1)*findpeaks(-1*h);
+        for i = 1:(length(hpeaks)-1)
+            if ((hpeaks(i) <= 0) || (hdips(i) >= 0) || (abs((hpeaks(i) + hdips(i))/2) >= thr))
+                temp = h;
+                test = 1;
+                break
+            end
+        end
+        k = k + 1;
+    end
+    y (n,:) = h;
+    % Step08
+    x = x - h;
+    n = n + 1;
+end
+
+end
+
+function drawgraph(y, t)
+subplot(3,1,1);				% Draw IMF01
+plot(t,y(1,:));
+title('IMF1');
+subplot(3,1,2);				% Draw IMF02
+plot(t,y(2,:));
+title('IMF2');
+subplot(3,1,3);				% Draw the Trend of Signal
+plot(t,y(3,:));
+title('Trend');
+end