1.0
Title Digital Filtering of Audio Signals
2.0
Objectives
To identify and analyze a digital audio signal with added noise.
To identify these unknown noise signal (filter specification) and to eliminate them using FIR as well as IIR digital filters (determining filter order and coefficient).
To apply the knowledge that you obtained from the lecture to accomplish the task by using MATLAB
3.0
Result
3.1
Part A: Comparison between Original and Noise Signal in Time Domain and Frequency Domain 1. Matlab Coding
Line
Command
Comment
1
[ori, fs]= wavread('D:\Benw\sem7\DSP\wave_files\G2o.
% get frequency sampling from original sound
wav'); 2 3 4
wavplay(ori,fs); fs; a0=0:1/fs:(length(ori)-1)/fs;
5 6
figure (1); subplot(2,2,1); plot(a0,ori);
7 8 9 10 11
title('Original signal in time domain'); xlabel('Time (s)'); ylabel('Amplitude'); grid; m=length(ori)-1;
12 13
t1=0:fs/m:fs; wavefft=abs(fft(ori));
14 15
figure(1); subplot(2,2,3); plot(t1,wavefft);
16 17 18 19 20
title('Original signal in frequency domain'); xlabel('Frequency (Hz)') ylabel('Magnitude'); grid; [noise, ft] = wavread('D:\Benw\sem7\DSP\wave_files\G2n.
% play the original sound from wavread % define the frequency sampling %takes depends on the original sound and the number of stages in the filter %create figure graphic object %define the subplot location (2,2,1) & plot the a0 and original signal %title of the subplot(2,2,1) %name of the label x as Time in second %name of the label y as Amplitude %put the grid on the graph %define the function n which is length of original song -1. %define the grid range depends on sampling size % fast fourier transform signal is equal absolute value %create figure graphic object %define the subplot location (2,2,3)& plot the t1 and fft signal %title of the plotted graph %x axis name as Frequency %y axis name as Magnitude %enable grid % to get frequency sampling from noise sound
wav'); 21 22 23
wavplay(noise,ft); ft; a0=0:1/ft:(length(noise)-1)/ft;
24
figure (1);
% play the noise sound from wavread % determine the frequency sampling %takes depends on the noise sound to be filtered and the number of stages in the filter %create figure graphic object
25
subplot(2,2,2); plot(a0,noise);
26 27 28 29 30
title('Noisy signal in time domain'); xlabel('Time (s)'); ylabel('Amplitude'); grid; n=length(noise)-1;
31 32
f2=0:ft/n:ft; wavefft2=abs(fft(noise));
33 34
figure(1); subplot(2,2,4); plot(f2,wavefft2)
35 36 37 38
title('Noisy signal in frequency domain'); xlabel('Frequency (Hz)'); ylabel('Magnitude'); grid;
%define the subplot location (2,2,2)& plot the a0 and noise signal %title of the plotted graph %x axis name as Time %y axis name as Magnitude %enable grid %define the function n which is length of noise song -1. %define the grid range depends on sampling size % fast fourier transform signal is equal absolute value %create figure graphic object %define the subplot location (2,2,4)& plot the f2 and fft signal %title of the plotted graph %x axis name as Frequency(Hz) %y axis name as Magnitude % allow the grid view
Figure 3.1.1: Comparison between Original Signal and Noise Signal in Time Domain and Frequency Domain.
Figure 3.1.2:Value of Sampling Frequency,ft.
Figure 3.1.3: The ways to find the Passband Frequency,fp and the Stopband Frequency, fs.
1. As the result in Matlab showed the frequency sampling,ft for both signal which are original and noise sound is 11025Hz. As in the Figure 3 below, we can see that the noise happen when 2240Hz, thus the chosen for passband frequency, fp is 2140Hz and the stopband frequency, fs is 2240Hz. So,
P
3.2
2f P 2 (2140 ) 0.388 ft 11025
The
cut-off
frequency,
fc
is
2190Hz.
C
( P S ) (0.388 0.406) 0.397 2 2
S
2f S 2 (2240 ) 0.406 ft 11025
As
f C ( f P f S ) / 2 2190 Hz
.
So,
Part B: FIR Filter Design and Filtering 1. Find the Order for FIR The requirement of the filter which the peak passband ripple is 3 dB and minimum stopband attenuation ois30 dB. Thus, the Hanning Window is more suitable choosing to design the filter, which in theoritical stopband attenuation for hanning window is αs ≤ 44dB. For the formula of Hanning are shows as below:
N Line 1
3.1 3.1 3.1 341 order f ( f S f P ) / f t (2240 2140 ) / 11025
Command [noise,freq]=wavread('D:\Benw\sem7\DSP\wave_fi
Comment % play the
noise sound
les\G2n.wav'); 2
a=fir1(341,0.397,hanning(342));
3 4 5
figure (4); subplot(3,2,1); stem (a);
6 7
grid; title('Impulse Response Coefficient');
% determine order = 341, cut-off frequency = 0.397 %create figure graphic object %define the subplot location % plot the discrete data of filter hanning window % allow the grid view %title of the plotted graph
8 9 10 11 12
xlabel ('Time Index n'); ylabel('Amplitude'); figure (4); subplot (3,2,3); [h,w] = freqz(a,1,512);
13 14 15 16 17
plot (w/pi, 20*log10 (abs(h))); grid; xlabel ('\omega/\pi'); ylabel ('Gain, dB'); title ('Lowpass filter designed Window'); axis ([0 1 -80 10]); y = filter(a,1,noise); wavwrite(y,'fir1.wav'); [signal,freq]=wavread('fir1.wav'); n = length(signal); sound (signal,freq);
18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
using
Hann
%name of the label x as Time Index n %name of the label y as Amplitude %create figure graphic object %define the subplot location %digital frequency response(numerator and denominator) % plotted to define range window % allow the grid view %name of the label x %name of the label y %title of the plotted graph %range of axis for Hanning window plot
% to write the data stored % to get frequency sampling from noise % define length of signal % play the noise filter file ‘fir1.wav’ from wavread figure (4); %create figure graphic object subplot (3,2,2); %define the subplot location plot((1:n)/freq,signal); % plotted filter signal in time domain title('Filtered signal in time domain'); %title of the plotted graph xlabel('Time (s)'); %name of the label x ylabel('Amplitude'); %name of the label y grid; % allow the grid view figure(4); %create figure graphic object signal0=signal-mean(signal); % the value of signal10 fsignal=fft(signal0); % the axis value of fsignal subplot (3,2,4); %define the subplot location plot((1:n/2)/n*freq,abs(fsignal(1:n/2))); % plotted filter signal in freq domain grid; % allow the grid view title('Filtered signal in frequency domain'); %title of the plotted graph xlabel('Frequency (Hz)'); %name of the label x ylabel('Magnitude'); %name of the label y figure(4); %create figure graphic object subplot (3,2,[5 6]); %define the subplot location zplane(a, 1) % plot zero pole at z-plane title('Z-plane'); %title of the plotted graph grid; % allow the grid view Table 1: Matlab Coding (design filter response coding)
Figure 3.2.1: Plot the designed filter response 341 order Line 1 2 3 4
Command [ori, fs] = wavread('D:\Benw\sem7\DSP\wave_files\G2o.wav' wavplay(ori,fs); fs; a0=0:1/fs:(length(ori)-1)/fs;
5 6 7 8 9 10 11
figure (1); subplot(3,2,1); plot(a0,ori); title('Original signal in time domain'); xlabel('Time (s)'); ylabel('Amplitude'); grid; m=length(ori)-1;
12 13
t1=0:fs/m:fs; wavefft=abs(fft(ori));
14 15
figure(1); subplot(3,2,2); plot(t1,wavefft)
16 17 18 19 20
title('Original signal in frequency domain'); xlabel('Frequency (Hz)') ylabel('Magnitude'); grid; [noise, ft] =
('D:\Benw\sem7\DSP\wave_files\G2n.wav'); 21 22 23
wavplay(noise,ft); ft; a0=0:1/ft:(length(noise)-1)/ft;
Comment % get frequency sampling from original sound % play the original sound from wavread % define the frequency sampling %takes depends on the original sound and the number of stages in the filter %create figure graphic object %title of the plotted graph %title of the subplot(2,2,1) %name of the label x as Time in second %name of the label y as Amplitude %put the grid on the graph %define the function n which is length of original song -1. % fast fourier transform signal is equal absolute value %create figure graphic object %define the subplot location (3,2,2)& plot the t1 and fft signal %title of the plotted graph %name of the label x %name of the label y % allow the grid view % to get frequency sampling from noise sound % play the noise sound from wavread % determine the frequency sampling %takes depends on the noise sound to be filtered and the number of stages in the filter
24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
figure (1); subplot(3,2,3); plot(a0,noise);
%create figure graphic object %define the subplot location (3,2,3)& plot the a0 and noise signal title('Noisy signal in time domain'); %title of the plotted graph xlabel('Time (s)'); %name of the label x as Time in second ylabel('Amplitude'); %name of the label y as Amplitude grid; % allow the grid view n=length(noise)-1; % length of the signal f2=0:ft/n:ft; % range of frequency f2 wavefft2=abs(fft(noise)); % fast fourier transform signal is equal absolute value figure(1); %create figure graphic object subplot(3,2,4); plot(f2,wavefft2) %define the subplot location (3,2,4)& plot the f2 and fft signal title('Noisy signal in frequency domain'); %title of the plotted graph xlabel('Frequency (Hz)') %name of the label x as Frequency in Hz ylabel('Magnitude'); %name of the label y as Magnitude grid; % allow the grid view n = length(signal); % define length of signal sound (signal,freq); % play the noise filter file ‘fir1.wav’ from wavread figure (1); %create figure graphic object subplot (3,2,5); %define the subplot location plot((1:n)/freq,signal); % plotted filter signal in time domain title('Filtered signal in time domain'); %title of the plotted graph xlabel('Time (s)'); %name of the label x as Time in second ylabel('Amplitude'); %name of the label y as Amplitude grid; % allow the grid view figure(1); %create figure graphic object signal0=signal-mean(signal); % the value of signal10 fsignal=fft(signal0); % the axis value of fsignal subplot (3,2,6); %define the subplot location plot((1:n/2)/n*freq,abs(fsignal(1:n/2))); %delay the truncated h[n] by n/2 samples grid; % allow the grid view title('Filtered signal in frequency domain'); %title of the plotted graph xlabel('Frequency (Hz)'); %name of the label x as Frequency in Hz ylabel('Magnitude'); %name of the label y as Magnitude Table 2: Comparison between Original Signal, Noise Signal and Filtered Signal
Figure 3.2.2:Comparison FIR design before and after filter noise.
Part C: IIR Filter Design and Filtering 1. Find the Order for IIR The requirement of the filter which require the peak passband ripple is 3 dB and minimum stopband attenuation is 30 dB. Thus, the Buterworth Window is more suitable choosing to design the filter. (1 P ) 2 1 d 2 S 1
1/ 2
(0.708) 2 1 2 (0.032) 1
1/ 2
0.995 975.563
1/ 2
0.032
P 20 log 10 (1 P ) dB, P 3dB (1 P ) 0.708 S 20 log 10 ( S ) dB, S 30dB ( S ) 0.032 2 0.388 P T 2 0.388 0.956 k 2 0.406 0.406 S T 2 log d log 0.032 N log k log 0.956
The order of the system integer required in estimation
N 70 order
Line
1
Command [noise,freq]=wavread('D:\Benw\sem7\DSP\wave_file
Comment % play the noise sound
s\G2n.wav' 2
[b,a]= butter(60,0.397);
3
[h, omega] =freqz(b,a,512);
% determine order= 10, cut-off frequency= 0.48 % digital frequency response(numerator and denominator)
4 5 6
y = filter(b,a,noise); wavwrite(y,'fir1IIR.wav'); [signal,freq]=wavread('fir1IIR.wav');
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
n = length(signal); sound (signal,freq); figure(6); subplot(3,2,1); plot(omega/pi, 20*log(abs(h))); grid; xlabel('\omega/\pi');ylabel('Gain,in dB'); ylabel('Gain,in dB'); title('Type Butter Filter'); axis([0 1 -60 5]); figure(6); subplot(3,2,3); plot(omega/pi, unwrap(angle(h))); grid; axis([0 1 -8 1]); xlabel('\omega/\pi'); ylabel('Phase,radians'); title('Type Butter Filter'); figure (6); subplot (3,2,2); plot((1:n)/freq,signal); title('Filtered signal in time domain'); xlabel('Time(s)'); ylabel('Amplitude'); grid; figure(6); signal0=signal-mean(signal); fsignal=fft(signal0); subplot (3,2,4); plot((1:n/2)/n*freq,abs(fsignal(1:n/2))); grid; title('Filtered signal in frequency domain'); xlabel('Frequency (Hz)'); ylabel('Magnitude'); figure(6) subplot (3,2,[5 6]); zplane(b, a); title('Z-plane'); grid; Table 3: Matlab Coding
% the parameters for filter y % to write the data stored % to get frequency sampling from noise sound % define length of signal % x and y of signal sound %create figure graphic object %define the subplot location % plotted filter signal in time domain % allow the grid view %name of the label x %name of the label y %title of the plotted graph % the scale of axis x and y %create figure graphic object %define the subplot location % to produce smoother phase plots % allow the grid view % the scale of axis x and y %name of the label x %name of the label y %title of the plotted graph %create figure graphic object %define the subplot location % the scale of axis x and y %title of the plotted graph %name of the label x %name of the label y % allow the grid view %create figure graphic object % the value of signal10 % the axis value of fsignal %define the subplot location % plotted filter signal in freq domain % allow the grid view %title of the plotted graph %name of the label x %name of the label y %create figure graphic object %define the subplot location % plot zero pole at z-plane %title of the plotted graph % allow the grid view
Figure 3.2.3: The Butterworth Designed Low Pass Filter 70th order Response
2.
Figure 3.2.4: Comparison IIR Design before and after filter noise.
4.0 DISCUSSION i. Plotting Audio Signal without noise and Audio Signal with Noise in Time and Frequency Domain Based on the filter knowledge, the MATLAB present a few built-in functions that allocate one to import and export audio files. It is understood that audio files sustained to be use in Matlab is Microsoft WAV ('pronounced 'wave') format. What we mean by built-in functions is, it is can be read and written to/from Matlab using built-in 'wavread' and wavwrite' functions.
To create the given audio signal in time and frequency domain, the command „wavread‟ is used followed by the location of the audio files, in bracket. If this line of command was skipped, there is no audio file will be retrieved. It is noted that the filenames and any text messages must be single quoted in Matlab. The command 'wavplay' is also require to be used in order to play the audio file. The frequency sampling from audio signal is obtained by using command 'fs'. The command [ori, fs] = wavread('audio file location') at line 1 will build the array ori contains the sound data and fs is the sampling frequency. The data is sampled at the same rate as that on a music CD (fs=44,100 samples/second). We say like this because, analog audio is recorded by sampling it 44,100 times per second and then these samples are used to recreate the audio signal when playing it back. To compute audio spectra in Matlab, the function 'fft' is used to computes FFT of that two audio signals. The magnitude/phase values of FFT coefficient was determined by using command 'abs' (line 13). The length of the plotted audio signal is written as t1=0:1/fs:(length(ori)-1)fs. From the result in Matlab showed the frequency sampling,ft for both original and noise sound is 11025Hz. As in the Figure 3, the passband and stopband edge frequency are identified by zooming the output signal which shows that fs=2240Hz and fp=2140Hz. Low pass filter and high pass filter were actually the symmetry of signal. Only half of the filter was taken because it is symmetry so low pass filter is chosen. Based on the fs and fp, the cut off frequency, ώc can be calculated. ii. FIR Filter Design In this task, the window technique is used for FIR filter design. The filter order and its coefficient are first determined. Given, peak passband ripple of 3 dB and minimum stopband attenuation of 30 dB. So, Hanning Window was choosen. The FIR order was 341 because according “a = fir(341,0.397,hanning (342))” command, which number of order is calculated due to N = 3.1/DF for Hanning window. While “hanning(342)” shows a causal FIR transfer function H(z) of length N+1 → 341+1= 342 is a polynomial in z-1 of degree N. While the 0.397=cut off frequency. Then, the function of
“[h,w] = freqz(b,1,512) “ command is to use the transfer
function
associated with the multirate filter to calculate the frequency response of the filter with the current coefficient values. The “wavwrite(y,Fs,N,filename)” command used to write the data stored in the variable y to a WAVE file called filename.
The command 'wavwrite(y,'fir1.wav'); is to save an array as a .wav audio file. This audio file is the new audio file where the noises have been filtered. It is automatically inside matlab folder path with
name fir1.wav. It depends on user, if they want to named it as 'mysong' so it should be written as : wavwrite(y,'mysong.wav').
To
get
frequency
sampling
from
noise
sound,
command
[signal,freq]=wavread('fir1.wav') was used. The “plot ((1:n/2)/n*freq, abs(fsignal(1:n/2)))” command especially n/2 used because h[n] is finite for positive and negative time for ideal filter is non-causal. To overcome the problem the solution is by delaying the truncated h[n] by n/2 samples. After FIR filter been applied (Figure 4), due to the z-plane in FIR filter design there are zeroes include in unit circle. In MATLAB
zplane(b, 1) command representing where b are the coefficients
obtained from fir1(),while a=1 since this is an FIR filter. The results based on Figure 4 showed the filtered signal graph, noise signal and original signal in the form of time domain and frequency domain. The signal obtained after filtering process by using Hanning window approximately similar with the original signal in time domain and frequency domain. However, in frequency domain the filtered signal sample has half of the sample compared to original signal‟s sample. This is because the resulting impulse response, h[n-M/2] is causal, stable FIR filter.
Based on comparison (Figure 5), there are obvious differences between audio file without noise and with noise, in time and frequency domain. In time domain, audio signal with noise looked almost periodic over short time interval compare to the audio signal without noise, which shows that original audio signal has very high frequency due to the presence of noise. While in frequency domain, the noises are clearly can be identified compared to time domain. In frequency domain, it shows that low frequency signal is being interrupted by high frequencies approximately at 2240Hz and the signal is symmetric. Therefore, the elimination of noise is implemented in frequency domain by removing high frequencies (noises that attached to original audio signal) to retrieve back the original audio signal. After the FIR filter been applied, a high frequency components have been successfully suppressed in the output, which satisfied our expectation and theory. The audio files seems to reduce its noisy sound when we playing the audio file.
iii. IIR Filter Design In this IIR filter design, given peak passband ripple of 3 dB and minimum stopband attenuation of 30 dB. So, Butterworth Window was choosen because Butterworth refers to a type of filter response. It is sometimes called the Maximally Flat approximation, because for a response of order n, the first
derivatives of the gain with respect to frequency are zero at frequency = 0. There is no ripple in the low pass band, and DC gain is maximally flat. This command representing [h, omega] =freqz(b,a,512), while freqz is used to declare b and a which contain numerator and denominator. The z-plane for IIR design filter showed that the poles and zeroes included in unit circle (Figure 6). In this IIR filter, the order,N is 60. The function of MATLAB [b,a]= butter(70,0.397) designs an th
N order=60 low pass digital Butterworth filter and returns the filter coefficient in length N+1 vectors b=numerator and a=denominator. The cut off frequency wn = Y.YY because wn must be 0.0< wn<1.0 with 1.0 corresponding to half the sample rate. The signal obtained after filtering process by using Butterworth filter is approximately similar compared to the original signal in the form of time domain and frequency domain.
Based on comparison (Figure 7), there are obvious differences between audio file without noise and with noise, in time and frequency domain. In time domain, audio signal with noise looked almost periodic over short time interval compare to the audio signal without noise, which shows that original audio signal has very high frequency due to the presence of noise. While in frequency domain, the noises are clearly can be identified compared to time domain. In frequency domain, it shows that low frequency signal is being interrupted by high frequencies approximately at 2240Hz and the signal is symmetric. Therefore, the elimination of noise is implemented in frequency domain by removing high frequencies (noises that attached to original audio signal) to retrieve back the original audio signal. After the IIR filter been applied, again, the same result occur as in FIR, whereby we able to notice the absence of the highest frequency in the output sample. iv. FIR vs. IIR FIR and IIR used to filter audio signal to removed or reduced noise. So, now we can summarize that which type of digital filter is the best. Based on our comparison, can obtain that the number of order for FIR is larger than IIR. (In calculation, Order FIR=342, Order IIR=60). Based on this, we can state that FIR filter need higher order and trouble-free to apply but FIR does not have feedback, which make FIR stable than IIR. Vice versa, IIR filter are difficult to implement but they need lower order and not stable. The reason IIR unstable is because IIR have poles in their transfer function which arise a tendency that the filter can become unstable. This can be verify in Figure 6, the pole-zero plot graph, we saw the poles may progress to the outside of the unit circle, and the zeros may spread around z=-1.
The FIR can be designed with exact linear phase, filter structure always stable with quantize coefficient because FIR have no poles but only have zeros. The purpose of phase linearity is to keep away from distortions in the output signal, which is in our case, it is an audio signals.
In review, IIR is infinite and used for application where linear characteristic are not of concern. While FIR are finite which required for linear-phase characteristics. FIR filters are often preferred over IIR because they are more stable and feedback is not involved.
5.0 CONCLUSION An audio file explains a format, sometimes referred to as the 'container format', for storing digital audio data. The audio signals in ".wav" format that was with noise, it was successfully filtered out by using filtering methods. There are two types of digital filtering techniques we are using, which are Finite Impulse Response (FIR) and Infinite Impulse Response (IIR). Both filter are suitable to reduce and removed noise. Firstly have to gather signal characteristics and well understood about the design process to perform FFT to know the frequency components in signals, design a filter where you can remove unwanted signals. After that apply filter to our signal to remove unwanted signal or noise. The noises are successfully removed. It retains the low frequency (bass) and soften the high frequency (treble). We say it as 'soften' because it is not exactly removed all the noises, it is only reduce the noises. As observation, the method of filtering noise is depending on the application of the problem.FIR is suitable for linear phase and stable output, While IIR are suits for always not linear but can be designed to be stable At the end of these assignments empower us to apply musical capacity of Matlab. We believe Matlab could be a decent digital music production tool. This is useful for musicians to further discover their music as well as engineer to better considerate music. To conclude, both FIR and IIR have its advantages and disadvantages. 6.0 REFERENCES
1. 2.
S.K. Mitra, Digital Signal Processing: A Computer-Based Approach, New York, NY: McGraw-Hill, 1998 http://en.wikipedia.org/wiki/Filter_(signal_processing)
3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13.
http://dsp.stackexchange.com/questions/9661/estimating-filtered-noise-variance Dr. DePiero, Filter Design by Frequency Sampling, CalPoly State University http://class.ee.iastate.edu/mmina/ee186/labs/Audio.htm W.James MacLean, FIR Filter Design Using Frequency Sampling A. V. Oppenheim. (2008). Signals and Systems. Prentice Hall. S. D. Stearns and R.A. David. Signal Processing Algorithms in MATLAB. Prentice Hall, 1996. Mathworks Inc., MATLAB Reference Guide, Mathworks, Natick, MA,2000 Lawrence R. Rabiner, Linear Program Design of Finite Impulse Response Digital Filters, IEEE 1972 Maurice G.Bellanger, Adaptive Digital Filters second edition, Marcel dekker 2001 http://www.cs.tut.fi/sgn/arg/intro/basics.html
