How do i filter sine wave which has noise of random numbers -0.5 to 0.5?

So this is the actual question.i need to select a filter and we are supposed to use the function filter (b,a,X) how does one define parameters for b and a? Calculate and plot a sine signal between 0 and 2.pi and a step size of 0.001, which is influenced with noise (random numbers between -0.5 and 0.5), plot and one subplot should include a filtered sine

thanks

Answer 1

[Geoff Hayes]

Rohin - is this a homework question?

There seems to be three parts to it - plot a sine signal in the interval [0,2pi] with a step size of 0.001; add some noise to the signal; and apply a filter to the noisy data and plot it.

Start with just the plotting of the sine signal and adding noise to it. See rand and in particular the example on how to generate random numbers in the interval [a,b] (which in your case, is [-0.5,0.5]).

Then look at filter and look at the one example. Try it out with your code, experimenting with different window sizes, plotting the results each time. Observe what happens as the window size increases.

Start with the above and see what happens!

[Rohin]

yes this is a homework question. thank you for your answer. i have done the previous questions and only asked help for the parts i cannot understand, the problem is i am not able to understand how the filter functions works and how a and b are determined or defined?. trial and error is giving me a sinewave but its still very noisy, and not a relatively smooth curve, however, trial and error doesnt help much if i cant understand the basic concept involved.

[Geoff Hayes]

Rohin - from the link for filter , the equation MATLAB uses to relate a, b, and x is as follows

 a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
                 - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

So let us look at their (easy) example where

 windowSize = 5;
 b = ones(1,windowSize)/windowSize;
 a = 1;

This means that the above equation reduces to

 y(n) = b(1)*x(n) + b(2)*x(n-1) + b(3)*x(n-3) + b(4)*x(n-4) + b(5)*x(n-5)
      = 1/5*x(n) + 1/5*x(n-1) + 1/5*x(n-3) + 1/5*x(n-4) + 1/5*x(n-5)
      = 1/5*(x(n) + x(n-1) + x(n-3) + x(n-4) + x(n-5))

as

 a(1) = 1, a(2) = 0, etc.

In this case, the output for any n is just the average of that value with the previous four elements (hence the description for this example as a running average). (Note that for n<5, it is not quite the average since less than five elements are used...)

If our a had more than one element in it, then the equation becomes a difference equation, making use of the outputs from previous iterations i.e. a(2)*y(n-1), a(3)*y(n-2), etc.

If you haven't gotten to the point in your class where you are designing filters, it may be that the intent of this question is to just realize what the effects of this filter are on a noisy sine wave. How does the window size (number of elements in b) affect the noisy input signal? How does a influence the output?

Answer 2

[Image Analyst]

I'd just filter with a Savitzky-Golay filter of order 4 or 5, which would be a pretty darn good match for a sine wave (if you remember the Taylor series form of the sine wave).

smoothY = sgolayfilt(y, polynomialOrder, windowWidth);

As you probably know, a Savitzky-Golay filter fits the data in a window to a polynomial, and then the window slides along your signal getting a fitted value for every location. In your question "a" and "b" would be the polynomialOrder, and the windowWidth.

Attached is a demo, which requires the Signal Processing Toolbox.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables. Or clearvars if you want.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 25;

% Make a noisy sine wave signal
x = 1 : 300;
period = 50
y = sin(2*pi*x/period);
noiseAmplitude = 0.8;
y = y + noiseAmplitude * rand(size(y));
subplot(2,1,1);
plot(x, y, 'b-', 'LineWidth', 2);
grid on;
title('Noisy Signal', 'FontSize', fontSize);
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off') 

% Now smooth with a Savitzky-Golay sliding polynomial filter
windowWidth = 27
polynomialOrder = 3
smoothY = sgolayfilt(y, polynomialOrder, windowWidth);
subplot(2,1,2);
plot(x, smoothY, 'b-', 'LineWidth', 2);
grid on;
title('Smoothed Signal', 'FontSize', fontSize);

[Rohin]

Hi thanks for your answer, ive seen this at another post as well :) however the problem is still the same. how do you determine the windowwidth is 27? and how do you determine the polynomial order is 3? same with filter(b,a,X) apart from X, how do i calculate a and b?

[Image Analyst]

As you know a sine wave is a1*x + a3*x^3 + a5*x^5 + .......

The higher the order fit you put in, the more it will follow your original signal, and consequently do less smoothing. You would not want to use an order of like 11 or 17 or something because it won't do any smoothing.

The smaller the window, the less elements it's taking into consideration to do the smoothing and so the more it will follow your original signal. The bigger the window, the more the average will be "true" rather than influenced by local noise. However if the window gets too big, like a half cycle or something, then the mean is not really accurate. If the window was more than a whole cycle, the mean would always be close to zero no matter where it was located. So your window should be some fraction of a period, like a 10th of it or so.

Use trial and error to get some parameters that work well.

Or you can do it with a band pass filter in the Fourier domain. Call fft, zero out anything not at that exact frequency, then call ifft().

Answer 3

[Star Strider]

If you need to use Y = filter(b,a,X), and you have the (significant) advantage of knowing the frequency of the signal you want to filter, the Signal Processing Toolbox can do everything for you. (As wonderful as sgolay is, it seems out of bounds in your assignment.)

You can certainly calculate a and b by writing your own code to do it, but assuming you have the Signal Processing Toolbox functions and you need to learn how to use them, I would:

1. Choose a bandpass filter (a Butterworth design would probably work well here) centred on the frequency of your signal, with appropriate bandwidth, and with a reasonably steep rolloff in the stopbands,
2. If you choose a Butterworth filter, use buttord to define the best filter order and passband,
3. Use those outputs and the appropriate filter design function (if you choose a Butterworth design, butter) to calculate the numerator and denominator coefficient vectors,
4. Check the design with the freqz function to be sure it’s stable. If it is, use the b and a coefficient vectors with filter to filter your signal. If it is not, go back to Step #2 and redesign it. (NOTE that the MATLAB filter design functions will suggest using second-order-section implementation of your filter. I would normally suggest this, but since you have to use the numerator and denominator coefficient vectors, you can’t use sos on this assignment.)

[Rohin]

thank you for your help. I was able to finally understand the proper use of the butterworth filter and then use filter(b,a,X). thanks

[Star Strider]

My pleasure!