Sunday, October 13, 2019

IIR Filter: Direct Form I with MATLAB GUI

IIR Filter: Direct Form I


Here, you will build a desktop application (UI, user interface) to simulate how is filtered using IIR filter with direct form 1. A number of discrete signals are used as test signals are impulse, step, real exponential, sinusoidal, random, square, angled triangle, equilateral triangle, and trapezoidal signals. In the UI, this IIR filter is used to filter signals and audio file (wav) as well.


Follow these steps below:


1. Type this command on MATLAB prompt:


guide

The GUIDE Quick Start dialog box will show up, as shown in figure below:



2. Select the Blank GUI (Default) template, and then click OK:


3. Drag all components needed, as shown in figure below:


4. Right click on DRAW button and choose View Callbacks and then choose Callback. Define pbDraw_Callback() function below to read signals parameters and displays chosen signal in axes1 component:


% --- Executes on button press in pbDraw.
function pbDraw_Callback(hObject, eventdata, handles)
% hObject    handle to pbDraw (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

global x;

% Reads n0, n1, and n2
n0 = str2num(get(handles.n0,'String'));
n1 = str2num(get(handles.n1,'String'));
n2 = str2num(get(handles.n2,'String'));
exp = str2num(get(handles.expCoef,'String'));

% Reads parameters for sinusoidal
A = str2num(get(handles.A,'String'));
Freq = str2num(get(handles.Freq,'String'));
Fase = str2num(get(handles.Phase,'String'));
Fs = str2num(get(handles.Fs,'String'));

% Reads parameters for square and triangle signals
width = str2num(get(handles.width, 'String'));

% Choose option from popup menu popupmenusignal1
switch get(handles.popupmenusignal1,'Value')  
    case 1
       n = [n1:n2]; x = [(n-n0) == 0];
       
       % Display discrete signal
       axes(handles.axes1)
       stem(n,x,'linewidth',2,'color','b'); title('Discrete Impulse Signal');

    case 2
       n = [n1:n2]; x = [(n-n0) >= 0];
       % Display discrete signal
       axes(handles.axes1)
       stem(n,x,'linewidth',2,'color','b'); title('Discrete Step Signal');
       
    case 3
       n = [n1:n2]; x = [(n-n0) >= 0].*[(exp).^(n-n0)];
       % Display discrete signal
       axes(handles.axes1)
       stem(n,x,'linewidth',2,'color','b'); 
       title('Discrete Real Exponential Signal');
       
    case 4
       n = [n1:n2]; x = [(n-n0) >= 0].*[A*sin(2*pi*(Freq)*((n-n0)/Fs)+Fase)];
       % Display discrete signal
       axes(handles.axes1)
       stem(n,x,'linewidth',2,'color','b'); 
       title('Discrete Sinusoidal Signal');
       
    case 5
       n = [n1:n2]; x = [(n-n0) >= 0].*(rand(1,(n2-n1+1))-0.5);
       % Display discrete signal
       axes(handles.axes1)
       stem(n,x,'linewidth',2,'color','b'); 
       title('Discrete Random Signal');

    case 6
       n = [n1:n2]; x = [(n-n0) >= 0];
       ny = [n1:n2]; xy = [(ny-width-n0-1) >= 0];
       x = x - xy;
       
       % Display discrete signal
       axes(handles.axes1)
       stem(n,x,'linewidth',2,'color','b'); title('Discrete Square Signal');

    case 7
       n = [n1:n2]; x = n.*[n >= 0];
       ny = [n1:n2]; xy = ny.*[(ny-width-1) >= 0];
       x = (x - xy)/width;
       nb = [n1+n0:n2+n0];
       
       % Reads signal range
       set(handles.n1,'string',(n1+n0));
       set(handles.n2,'string',(n2+n0));

       % Display discrete signal
       axes(handles.axes1)
       stem(nb,x,'linewidth',2,'color','b'); 
       title('Discrete Angled Triangle Signal');
 
    case 8
       n = [n1:n2]; x = n.*[n >= 0];
       x2 = [zeros(1,width), x(1:end-width)];
       x1 = -x;
       x1 = [zeros(1,0.5*width), x1(1:end-0.5*width)];
       nb1 = n1+n0; nb2 = n2+n0;
       nb = [nb1:nb2];
       
       % Read signal range
       set(handles.n1,'string',(nb1));
       set(handles.n2,'string',(nb2));

       % Returns discrete equilateral triangle signal
       x = (x + 2*x1+x2)/(0.5*width);
       
       % Display discrete signal
       axes(handles.axes1)
       stem(nb,x,'linewidth',2,'color','b'); 
       title('Discrete Equilateral Triangle Signal');
       
    case 9
       n = [n1:n2]; x = n.*[n >= 0];
       x2 = [zeros(1,width), x(1:end-width)];
       x3 = [zeros(1,3*width), x(1:end-3*width)];
       x1 = x;
       x1 = [zeros(1,2*width), x1(1:end-2*width)];
       nb1 = n1+n0; nb2 = n2+n0;
       nb = [nb1:nb2];
       
       % Reads signal range
       set(handles.n1,'string',(nb1));
       set(handles.n2,'string',(nb2));

       % Returns discrete trapezoidal signal
       x = (x - x2 - x1 + x3)/width;
       
       % Display discrete signal
       axes(handles.axes1)
       stem(nb,x,'linewidth',2,'color','b'); 
       title('Discrete Trapezoidal Signal');
end  

5. You can run the UI. Choose one of signals and then click DRAW button:


6. Right click on FREQUENCY RESPONSE OF FILTER button and choose View Callbacks and then choose Callback. Define pbFreqResponse_Callback() function below to calculate FFT of designed IIR filter:


% --- Executes on button press in pbFreqResponse.
function pbFreqResponse_Callback(hObject, eventdata, handles)
% hObject    handle to pbFreqResponse (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

global x;

% Reads b0 to b4
b0 = str2num(get(handles.b0,'String'));
b1 = str2num(get(handles.b1,'String'));
b2 = str2num(get(handles.b2,'String'));
b3 = str2num(get(handles.b3,'String'));
b4 = str2num(get(handles.b4,'String'));

% Reads a1 to a4
a1 = str2num(get(handles.a1,'String'));
a2 = str2num(get(handles.a2,'String'));
a3 = str2num(get(handles.a3,'String'));
a4 = str2num(get(handles.a4,'String'));

% Filter coeffs
b = [b0 b1 b2 b3 b4];
a = [1 a1 a2 a3 a4];

% Menghitung Tanggapan Frekuensi
[h,w] = freqz(b,a,'whole',2001);

axes(handles.axes2); 
plot(w/pi,20*log10(abs(h)),'linewidth',2,'color','r'); grid on;

ax = gca;
ax.YLim = [-100 20];
ax.XTick = 0:.5:2;
xlabel('Normalized Frequency  (\times\pi rad/sample)')
ylabel('Magnitude (dB)'); title('Frequency Response of Designed Filter')

7. You can run the UI. Choose one of signals and click on DRAW button. Then click on FREQUENCY RESPONSE OF FILTER button to calculate FFT of designed IIR filter:


8. Right click on FILTERING SIGNAL button and choose View Callbacks and then choose Callback. Define pbFilteringSignal_Callback() function below to read input signal and calculates filtered signal in axes3 component:


% --- Executes on button press in pbFilteringSignal.
function pbFilteringSignal_Callback(hObject, eventdata, handles)
% hObject    handle to pbFilteringSignal (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

global x;

% Reads b0 to b4
b0 = str2num(get(handles.b0,'String'));
b1 = str2num(get(handles.b1,'String'));
b2 = str2num(get(handles.b2,'String'));
b3 = str2num(get(handles.b3,'String'));
b4 = str2num(get(handles.b4,'String'));

% Reads a1 to a4
a1 = str2num(get(handles.a1,'String'));
a2 = str2num(get(handles.a2,'String'));
a3 = str2num(get(handles.a3,'String'));
a4 = str2num(get(handles.a4,'String'));

% Filter coefficients
b = [b0 b1 b2 b3 b4];
a = [1 a1 a2 a3 a4];

% Calculates filtered signal
y = filter(b,a,x);

axes(handles.axes3); 
t = 0:length(y)-1;  

stem(t,y,'linewidth',3,'color','g'); title('Filtered Signal')



9. Right click on LOAD AUDIO FILE button and choose View Callbacks and then choose Callback. Define pbLoadAudio_Callback() function below to read audio signal and display it in axes2 component:


% --- Executes on button press in pbLoadAudio.
function pbLoadAudio_Callback(hObject, eventdata, handles)
% hObject    handle to pbLoadAudio (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

global x_audio;

[filename,path] = uigetfile({'*.wav'},'load audio file');
[x,Fs] = audioread([path '/' filename]);
handles.x = x ./ max(abs(x));
handles.Fs = Fs;
axes(handles.axes2);
time = 0:1/Fs:(length(handles.x)-1)/Fs;
plot(time, handles.x);
axis([0 max(time) -1 1]); title('Audio Signal');

axes(handles.axes3);
specgram(handles.x, 1024, handles.Fs);
title('Signal Spectrum');
handles.fileLoad = 1;

handles.fileNoise = 0;
handles.fileFinal = 0;

set(handles.freqSamplingAudio, 'String', num2str(Fs));
set(handles.NSamples, 'String', num2str(length(handles.x)));

x_audio = handles.x;
guidata(hObject, handles);

10. You can run the UI. Choose one of audio file by clicking LOAD AUDIO FILE button. The audio signal will be displayed in axes2 component and signal spectrum will be displayed in axes3 component:


11. Right click on FILTERING AUDIO FILE button and choose View Callbacks and then choose Callback. Define pbFilteringAudio_Callback() function below to read audio signal and display it in axes2 component:


% --- Executes on button press in pbFilteringAudio.
function pbFilteringAudio_Callback(hObject, eventdata, handles)
% hObject    handle to pbFilteringAudio (see GCBO)
% eventdata  reserved - to be defined in a future version of MATLAB
% handles    structure with handles and user data (see GUIDATA)

global x_audio;

% Reads b0 to b4
b0 = str2num(get(handles.b0,'String'));
b1 = str2num(get(handles.b1,'String'));
b2 = str2num(get(handles.b2,'String'));
b3 = str2num(get(handles.b3,'String'));
b4 = str2num(get(handles.b4,'String'));

% Reads a1 to a4
a1 = str2num(get(handles.a1,'String'));
a2 = str2num(get(handles.a2,'String'));
a3 = str2num(get(handles.a3,'String'));
a4 = str2num(get(handles.a4,'String'));

% Filter coefficients
b = [b0 b1 b2 b3 b4];
a = [1 a1 a2 a3 a4];

% Calculates filtered audio
y = filter(b,a,x_audio);

axes(handles.axes2); 
t = 0:length(y)-1;  
plot(t,y);title('Filtered Audio')

axes(handles.axes3)
specgram(y, 1024, handles.Fs);
title('Filtered Signal Spectrum');
handles.fileDimuat = 1;

handles.x=y;
guidata(hObject, handles);

12. Now, you can run the UI and click FILTERING AUDIO FILE button and see the filtered audio and its spectrum:





No comments: