I have been writing a very simple code to remove noise from a signal. The signal is just a sinusoidal wave, the noise is a random matrix, and the noisy signal is the addition of both.
The code is:
close all;clear;clc;
%% Declarations
ts = 0.001;
fs = 1/ts;
fc = 5;
t = 0:ts:2;
Wn = pi*fc/(2*fs);
n = 3;
%% Preparation
signal = cos(2*pi*fc*t);
noise = rand(1, length(signal)); % Generating Random Noise
noisySignal = signal + noise;
%% Filtering Stage
[b,a] = butter(n, Wn, 'low');
filteredSignal = filter(b, a, noisySignal);
filteredSignal = filteredSignal - mean(filteredSignal); % Subtracting the mean to block DC Component
%% Plotting
figure(1)
subplot(3,1,1)
plot(t, signal, 'linewidth', 1.5)
title('Signal')
ylim([-1.5 1.5])
grid minor
subplot(3,1,2)
plot(t, noise)
title('Noise')
ylim([-1.5 2])
grid minor
subplot(3,1,3)
plot(t, noisySignal)
title('Noisy Signal')
ylim([-1.5 1.5])
grid minor
figure(2)
plot(t, filteredSignal, 'r', 'linewidth', 1.5)
hold on
plot(t, signal, 'linewidth', 1.5)
hold off
legend('Filtered Signal', 'Original Signal')
grid minor
ylim([-1.5 1.5])
Figure 2; which is the figure for comparing both the filtered signal and the original signal; always appear like the image below.
I believe that the Wn
variable is not right, but I don't know how to calculate the correct normalized frequency.
Following this example form Matlab's documentation, if you want the cutoff frequency to be at fc
Hz at a sampling frequency of fs
Hz, you should use:
Wn = fc/(fs/2);
[b,a] = butter(n, Wn, 'low');
However you should note that this will produce a Butterworth filter with an attenuation of 3dB at the cutoff frequency. Since your sinusoidal signal is generated at a frequency fc
, the filtered sinusoidal would have an amplitude of roughly 70% of the original signal:
If you want less signal attenuation, you should increase the filter cutoff frequency. Of course doing so will also let a bit more noise through, so the exact amount is a trade-off between on how much signal attenuation your application can tolerate and how much noise you need to get rid of. For example, adding a margin of 1Hz and increasing the filter order (which gives you less attenuation for the same margin) with
Wn = (fc+1)/(fs/2);
n = 7;
[b,a] = butter(n, Wn, 'low');
would give you: