Signal Processing for Beginners using MATLAB

We will cover the following topics:

• Signal Generation and Spectrum Analysis;
• Signal Filtering using FIR and IIR Filters;
• Smoothing Signal.

Signal processing is the analysis of signals that focuses on removing noise, upsampling, downsampling, extracting information, smoothing, filtering, synthesising or detect features of interest in a measured signal.

Why is signal processing required?

As we see, sensors are used everywhere: in mobiles, health monitoring devices, vehicles, IoT, communications systems etc.

Signals need to be processed so that the information that they contain can be extracted, analyzed, or converted to another type of a signal so we can use it. Signals can be temperature, pressure, humid, digital data. Any kind of sensor can produce signals.

Signal generation and visualization in time domain

Let’s assign sampling frequency 1500 Hz, duration of signal 0 to .1 seconds and signal frequency is .60.

Fs = 1500;                  % Sampling frequency 1500 Hz 
t = 0:1/Fs:.1;               % 0 to .1 sec duration
f = 60;                     % Frequency of signal 60 Hz


% Square Wave
squareWave = square(2*pi*f*t);
subplot(2,2,1)
plot(squareWave);
axis tight
title('Square Wave')

% Sawtooth 
sawtoothWave = sawtooth(2*pi*f*t); 
subplot(2,2,2)
plot(sawtoothWave);
axis tight
title('Sawtooth')

% Sine waveform 
sineWave = sin(2*pi*f*t); 
subplot(2,2,3)
plot(sineWave); 
axis tight
title('Sine Wave')

% Chirp Signal 
f0=0; % Instantaneous frequency at time 0 is f0
t1=.1; % At the time instances defined in array t
f1=200; % Instantaneous frequency at time t1 is f1 in Hz.
chirpSignal = chirp(t, f0,t1,f1);
subplot(2,2,4)
plot(chirpSignal);
axis tight
title('Chirp Signal')

Signal generation and visualization in frequency domain

We will plot a spectrogram of signals. A spectrogram tells us how frequency is varying with the respect to time. A spectrogram is also useful to identify if any noise is introduced with the respect to time.

• Divide the signal into sections of length 120, windowed with a Hamming window.
• Specify 110 samples of overlap between adjoining sections.
• Evaluate the spectrum at 120/2+1=61 frequencies and (length(t)−110)/(120−110)=139 time bins.

clear,clc % Clear variables and command window
Fs = 1500;                  % Sampling frequency 1500 Hz
t = 0:1/Fs:1;               % 0 to 1 second time duration
f = 60;                     % 60 Hz signal 

%Square Wave
squareWave = square(2*pi*f*t);
subplot(2,2,1)
spectrogram(squareWave,120,110,120,'yaxis')
title('Square Wave' )

%Sawtooth wave 
sawtoothWave = sawtooth(2*pi*f*t); 
subplot(2,2,2)
spectrogram(sawtoothWave,120,110,120,'yaxis')
title('Sawtooth')

%Sine waveform 
sineWave = sin(2*pi*f*t); 
subplot(2,2,3)
spectrogram(sineWave,120,110,120,'yaxis')
axis tight
title('Sine Wave')

f0=0; % Instantaneous frequency at time 0 is f0
t1=.1; % At the time instances defined in array t
f1=200; % Instantaneous frequency at time t1 is f1 in Hz.
chirpSignal = chirp(t, f0,t1,f1);
subplot(2,2,4)
spectrogram(chirpSignal,120,110,120,'yaxis')

axis tight
title('Chirp Signal')

Filter designing and signal filtering

Some points to be remembered while designing filters:

• When designing a filter, the first choice you make is whether to design an FIR or IIR filter.
• Choose FIR filters when a linear phase response is important.
• If the desired frequency response can be approximated by a smaller number of coefficients, then we choose IIR filters, where the phase linearity is not a concern.
• IIR filters consist of zeros and poles, and require less memory than FIR filters, whereas FIR filters consist only of zeros. So they require high memory because the order of FIR filters increases.
• Filter Designs – Specifying the Filter Order (Two Scenario)
  – If you are targeting a hardware which has constrained the filter order to a specific number.
  – Another common scenario is when you have computed the available computational budget (MIPS) for your implementation and this affords you a limited filter order.

Design FIR and IIR filters for a noisy signal where signal frequency is 100 Hz

Objective
– Compare IIR and FIR filters for the same noisy signal.
– Find the cost of both filters.

Create a signal and add noise

% Signal
clear,clc
Fs = 1000; % Sampling frequency 
t = 0:1/Fs:2; % Sample time from 0 to 2 seconds
x = cos(2*pi*100*t); % Signal
plot(t(1:100),x(1:100)) % Plotting 100 samples only

Add normal random noise to the original signal

noise=0.5*randn(size(t)) % Normal random noise
x=x+noise;
plot(t(1:100),x(1:100))

Visualize a noisy signal in frequency domain

We will use pwelch to convert a time domain signal to the frequency domain. It estimates welch’s power spectral density.

[pxx,fxx] = pwelch(x,window,noverlap,f,fs)

• Output pxx is power at frequency fxx.
• x is noisy signal.
• window uses the input vector or integer, window, to divide the signal into segments.
• noverlap uses noverlap samples of overlap from segment to segment.
• f returns the two-sided Welch PSD estimates at the frequencies specified in the vector, f.
• fs is sampling frequency.

To extract peaks we are going to use findpeaks, it extracts peaks from the signal, you can assign any number of peaks you want, in my case I’ve assigned a threshold value of height of the peak.

[pxx,fx] = pwelch(x,[],[],[],Fs); % For frequency domain
plot(fx,10*log10(abs(pxx))) % power is transformed in DB      
[peak, loc]=findpeaks(pxx,fx,'MinPeakHeight',.1) % Extraxt signal peak
hold on
plot(loc,10*log10(peak),'ro')
title(['Siginal frequency is ', ' ',num2str(loc)]) % title num2str convert numeric value in string.
hold off

As it is shown in the figure above, the signal frequency is almost 100Hz. We are going to design a low pass filter with cuttoff frequency of 120 Hz.

Design a FIR low pass filter

Design the filter using designfilt. Set the filter response to ‘lowpassfir‘ and input the specifications as Name, Value pairs. With designfilt, you can specify your filter design in Hz. Design a lowpass FIR filter for data sampled at 1000 Hz. The passband-edge frequency is 120Hz. Constrain the filter order to 20.

We are going to use a window method to design a low pass filter and the window type is kaiser with beta value 3.

Fs = 1000;
Hd = designfilt('lowpassfir','FilterOrder',20,'CutoffFrequency',120, ...
'DesignMethod','window','Window',{@kaiser,3},'SampleRate',Fs);

To visualize the filter magnitude response we will use fvtool.

fvtool(Hd)

Apply the filter to the signal and plot the result for the first ten periods of the 100 Hz sinusoid.

y = filter(Hd,x) filters the input data x using Hd filter.

y = filter(Hd,x);
plot(t(1:100),x(1:100))
hold on
plot(t(1:100),y(1:100))
xlabel('Time (s)')
ylabel('Amplitude')
legend('Original Signal','Filtered Data')

Using DSP Filter

The filter which we designed above is based on the signal processing toolbox in MATLAB. There is another way to design a digital filter if we know the coefficients of the filter.

To extract the coefficients of the filter we use tf function and to design a dsp filter we use dsp.FIRFilter

NUM = tf(Hd);

LP_FIR = dsp.FIRFilter('Numerator',NUM);

To apply the filter to data, we can use the filter command or we can directly input the noisy signal to the filter. To visualize the spectrum we use dsp.SpectrumAnalyzer.

LP_FIR = dsp.FIRFilter('Numerator',NUM); % Design dsp filter
SA     = dsp.SpectrumAnalyzer('SampleRate',Fs,'SpectralAverages',5); % System object for spectrum 

    y = LP_FIR(x'); % apply filter on noisy data
    step(SA,y);  % display spectrum

Plot the result for the first ten periods of the 100 Hz sinusoid.

figure
    plot(t(1:100),x(1:100))
hold on
plot(t(1:100),y(1:100))
xlabel('Time (s)')
ylabel('Amplitude')
legend('Original Signal','Filtered Data')

Filtered signal in frequency domain

[pxx,fx] = pwelch(y,[],[],[],Fs);
figure,
plot(fx,20*log10(abs(pxx)))
hold on
[pxx,fx] = pwelch(x,[],[],[],Fs);

plot(fx,20*log10(abs(pxx)))
hold off
legend('Filtered signal','Noisy Signal')

Designing IIR Filters for the same noisy signal

The filter order you obtain for an IIR filter is much smaller than the order of the corresponding FIR filter. That’s why IIR filters require low memory. In this case we are going to use filter order 2 only. In the FIR filter it was 20.

lpFilt = designfilt('lowpassiir','FilterOrder',2, ...
         'PassbandFrequency',120,...
         'SampleRate',1000); % Design a filter using signal processing
fvtool(lpFilt)            % Filter magnitude response
   NUMIIR = tf(lpFilt); % coefficients of IIR FIlter
   LP_IIR = dsp.FIRFilter('Numerator',NUMIIR); % Design a digital filter

Apply the IIR filter to the signal and plot the result for the first ten periods of the 100 Hz sinusoid.

SA = dsp.SpectrumAnalyzer('SampleRate',Fs,'SpectralAverages',5);
 y = LP_IIR(x'); 
    step(SA,y);   

Plot the result for the first ten periods of the 100 Hz sinusoid.

plot(t(1:100),x(1:100))
hold on
plot(t(1:100),y(1:100))
xlabel('Time (s)')
ylabel('Amplitude')
legend('Original Signal','Filtered Data')

Compare the FIR and IIR designs

fvtool(LP_FIR,LP_IIR,'Fs',Fs);

Compute the cost of both FIR and IIR filters

cost_FIR = cost(LP_FIR)
cost_IIR = cost(LP_IIR)

Smoothing Signal

There are different ways you can smooth your signal. They are given below. We are going to use moving mean method.

Plot original signal

We’re going to smooth solar intensity data.

plot(intensity)

The smoothed data using the moving mean method with smoothing factor 0.25.

% Smooth input data
smoothedData = smoothdata(intensity,'movmean','SmoothingFactor',0.25);

% Display results
clf
plot(intensity,'DisplayName','Input data')
hold on
plot(smoothedData,'r','LineWidth',1.5,...
    'DisplayName','Smoothed data')
hold off
legend

If you don’t want to write codes for data pre-processing then you can use the capability of live script.

If you want to see the capability of live script you can click here for the recorded webinar.

Or click here for the recorded webinar on Signal Processing with MATLAB.