Chiusakpung's Electronics Blog

%% = Bandpass Filter:
clear, clc
freq = [1, 2, 3, 4, 5];
angfreq = 2 * pi * freq;
t = 0 : 0.001 : 10;
x = sin(angfreq(1) * t) + sin(angfreq(2) * t) + sin(angfreq(3) * t) + sin(angfreq(4) * t) + sin(angfreq(5) * t);
[num, den] = butter(3, [0.1, 0.2]);
y = filter(num, den, x);
subplot(2,1,1); plot(t,x);
subplot(2,1,2); plot(t,y);

%% DFT using FFT function

clear, clc
t = 0:0.001:0.6;
x = sin(2*pi*50*t)+sin(2*pi*120*t);
y = x + 2*randn(size(t));
subplot(2,1,1); plot(1000*t(1:50),y(1:50));
title('Signal Corrupted with Zero-Mean Random Noise');
xlabel('time (milliseconds)');

Y = fft(y,512);
Pyy = Y.* conj(Y) / 512;
f = 1000*(0:256)/512;
subplot(2,1,2); plot(f,Pyy(1:257));
title('Frequency content of y');
xlabel('frequency (Hz)');

%% DFT using FFT function - sin & cos
% Samping Rate
% Sampling frequency fs
% > twice highest frequency of input
% signal! )
clear, clc
fs = 100; Ts = 1/fs; Nf = 1000;
t = 0:Ts:10;
x = sin(2 * pi * 2 * t);
X = fft(x,Nf);
Pxx = X.*conj(X);
f = fs * (0:(Nf/2))/Nf;
subplot(4,1,1); plot(t,x); title('Input signal x(t)'); xlabel('time t(s)');
subplot(4,1,2); plot(f, Pxx(1:(Nf/2) + 1)); title('Output signal X(f) - Fourier transform of x(t)'); xlabel('freq f(Hz)');
x2 = ifft(X, Nf);
subplot(4,1,3); plot((1:Nf)*Ts,x2); title('x_2(t) - Inverse Fourier Transfrom of X(f)');
subplot(4,1,4); plot((1:Nf)*Ts, x(1:Nf) - x2); axis([1*Ts, Nf*Ts, -1,1]); title ('Difference betw. x(t) & x_2(t)');

%% Creation of Transfer Function of BPF by tf & bode
clear, clc
wL = 2 * pi * 1000;
wH = 2 * pi * 5000;
gain = 10;
H = tf([1 0], [1/251340, 1/gain, 785.4003]);

%subplot(2,2,1); plot(f, H); % plot real(H)
%subplot(2,2,2); plot(f, H .* conj(H));
%subplot(2,2,3); plot(f, abs(H)); % plot |H(f)|
%subplot(2,2,4); plot(f, angle(H)); % plot ∠H(f)
%subplot(2,1,2);
bodeH = bodeplot(H);
setoptions(bodeH,'FreqUnits','Hz');
grid;
bandwidth(H)

%% Bandpass Filter
clear, clc
% Variables for Bode plotting - Range of freqeuncy (Hz) & Gain
Lorder = 1;
Horder = 5; 
gain = 10;

% Domain Setting
tmax = 0.05; tstep = 0.00001;
fmax = 50000; fstep = 1;
t = - tmax : tstep : tmax;
f =  (10^Lorder) - 9 : (10^Horder);  % - 9 for even number points
w = 2 * pi * f;

 % Nf  = Number of points for FFT         = size of X(f)  =  size of
 % frequency domain 'f' = size of H(f)
Nf = length(f);   
  

% for Plotting X(f) if dScaling = 3/5, the domain of frequency  : 1 ~ 50000
% * 3/5 = 1~30000 Hz
dScaling = 3/5;

% Transfer Function
% Gain : 10
% Cut-off 3dB freq : low = 1 kHz, high = 5 kHz
wL = 2 * pi * 1000;
wH = 2 * pi * 5000;

H_tf = tf([1 0], [1/251340, 1/10, 785.4003]);
H = (j * w) ./ (-w.*w/251340 + j*w/10 + 785.4003);

figure(1);
bodeH = bodeplot(H_tf);
setoptions(bodeH,'FreqUnits','Hz'); grid;

figure(2);
subplot(2,1,1); semilogx(w / (2 * pi), 20*log10(abs(H))); axis([10^Lorder, 10^Horder, -60, 20]); grid;  % plot Mag.
title('Bode plot of H(jw) - Not created by tf function');
ylabel('Magnitude (dB)');
subplot(2,1,2); semilogx(w / (2 * pi), unwrap(angle(H)) * 180 / pi); axis([10^Lorder, 10^Horder, -90,90]); grid; % plot Angle
ylabel('Angle (deg)');
xlabel('Freqeuncy (Hz)');
set(gca, 'YTick', -90 : 45 : 90);
set(gca, 'YTickLabel',{'-90','-45','0','45','90'});

% Initializing test signal x(t)
x = sin(2*pi*100*t) + sin(2*pi*200*t) + sin(2 * pi * 500*t);
x = x + sin(2*pi*1000*t) + sin(2*pi*2000*t) + sin(2*pi*5000*t);
x = x + sin(2*pi*10000*t ) + sin (2*pi*20000*t);

% X(f) = FFT {x(t)}s
X = fft(x, Nf);

% Output signal Y(f) = H(f)X(f) <-> y(t) = h(t) * x(t)
%size(H), size(X)
Y = H .* X;

MagX = abs(X);
MagY = abs(Y);

figure(3);
subplot(2,2,1); plot(t, x); title('Test Signal x(t)'); xlabel('time t (s)');
subplot(2,2,2); plot(w/(2*pi), abs(H)); title('Transfer Function H(f)'); xlabel('freqeuncy f (Hz)'); ylabel('Mag. (Not dB)');
subplot(2,2,3); plot(f(1: Nf/2 * dScaling), MagX(1: Nf/2 * dScaling)); title('Not filtered signal |X(f)| - Fourier Transformed'); xlabel('frequency f (Hz)'); grid;
SeeFreq_100_Hz = 0; % 1 for yes, 0 for no
subplot(2,2,4); plot(f(1 : Nf/2 * dScaling), MagY(1 : Nf/2 * dScaling)); title('Fitered signal |Y(f)| = |H(f)||X(f)|'); xlabel('frequency f (Hz)'); grid;
xlim([1,25000*(1 - SeeFreq_100_Hz) + 2500 * SeeFreq_100_Hz]);
set(gca, 'YTick', [50000/11, 50000/6, 50000/5, 50000/sqrt(7), 35355, 50000]);
set(gca, 'YTickLabel',{'50000/11', '50000/6', '50000/5', '50000/√7','50000/√2','50000 (Gain = x10)'});

% Inverse Fast Fourier Transforms
y = ifft(ifftshift(Y));
tScaling = 1/Nf;

figure(4);
Magy = real(y);
subplot(2,1,1); plot(tScaling * (0:length(Magy) - 1), Magy); title('Output Signal y(t)'); xlabel('time t(s)');axis([0,0.005,-25,25]);

xfilt = 50000/11 * sin(2*pi*100*t) + 50000/6*sin(2*pi*200*t) + 50000/sqrt(7)*sin(2 * pi * 500*t);
xfilt = xfilt + 35355 * sin(2*pi*1000*t) + 50000*sin(2*pi*2000*t) + 35355 * sin(2*pi*5000*t);
xfilt = xfilt + 50000/sqrt(7) *sin(2*pi*10000*t ) + 50000/5*sin (2*pi*20000*t);
xfilt = xfilt/50000;
xfilt = xfilt * 10; % Gain = 10;

subplot(2,1,2); plot(t,xfilt); title('Expected Output Signal y(t)'); xlabel('time t(s)'); axis([0,0.005,-25,25]);

%close(1);
%close(2);
%close(3);
%close(4);