DSP: Using a basic digital filter: Stock market example

Digital filters are very powerful tools in DSP. In this example we will use a digital filter to help with stock market analysis. In the figure above, you can see the raw data plotted with a 5 day moving average and a 20 day moving average. These two moving average lines can be used for buy/sell predictions.

To calculate the 20 day moving average, the window size and filter coefficients must be established. The following Matlab code sets the window size and calculates the filter results. The five day filter is implemented using the same code, with a windowSize = 5. The matrix x in the example contains the raw stock data. For more info on the filt function: http://www.mathworks.com/access/helpdesk/help/techdoc/ref/filter.html

windowSize = 20; % Twenty day average
% The following filter applies the 1/windowSize coefficients
twentyDay = filter(ones(1,windowSize)/windowSize,1,x);

Use of Windowing Techniques on 3G CDMA Signal

Here is an example of how windowing (blackman in this example) can help reduce sidebands. This technique is useful in making spectral measurements such as Adjacent Channel Leakage Ratio. By minimizing the sidebands generated by an arbitrary waveform generator, for example, the true device performance can be more accurately measured. For more information on windowing, see the following article from National Instruments: http://zone.ni.com/devzone/cda/tut/p/id/4844

Coherent Noise Averaging

% James Eastham
% Member, IEEE
% Revision: R1
% Contact: james.eastham@ieee.org
%
% Description: This script shows the affect of coherent averaging to
% reduce random noise.

%% Affect of Additive Noise
close all; % close all plots
clear all; % clear all variables
fs = 500; % sample freq
Ts=1/fs; % sample period
n=0:1*fs; % time window (fs=samples per second), 1 second
tn=n*Ts; % sample index
x=cos(2*pi*50*tn); % 50Hz sinusoid
k=.5; %noise amplitude
n=randn(size(x)); % random noise
s=x+(k*n); %signal + random noise with amplitude k
% Time domain plot of x(t)=cos(2*pi*50*tn)
figure('Color',[1 1 1]);
plot(tn,s,'g','LineWidth',2);box off;hold on;
title('Sinusoid w/ Additive White Noise');
xlabel('Time(s)');
ylabel('Amplitude');
% Frequency Domain Plot
figure('Color',[1 1 1]);
[X,f]=ComputeSpectrum(s,fs,2^12);
orig_X=X;
h=plot(f,20*log(X));
box off;
grid on;
set(h,'Linewidth',1);
xlabel('Frequency(Hz)');
ylabel('Amplitude(dB)');
title(['Frequency Spectrum w/ Additive Noise k=.5']);
axis([0 250 -50 120]);

%% Example of Noise Averaging, 3 Coherent Averages
fs = 500; % sample freq
Ts=1/fs; % sample period
n=0:1*fs; % time window (fs=samples per second), 1 second
tn=n*Ts; % sample index
x=cos(2*pi*50*tn);
k=.2; %Noise amplitude
n1=randn(size(x));
n2=randn(size(x));
n3=randn(size(x));
s1=x+(k*n1);
s2=x+(k*n2);
s3=x+(k*n3);
average_s=((s1+s2+s3)/3);
% Time domain plot of x(t)=cos(2*pi*50*tn)
figure('Color',[1 1 1]);
h=plot(tn,average_s,'g','LineWidth',2);box off;hold on;
title('Sinusoid w/ Additive White Noise: 3 Averages');
xlabel('Time(s)');
ylabel('Amplitude');
% Frequency Domain Plot
figure('Color',[1 1 1]);
[X,f]=ComputeSpectrum(average_s,fs,2^12);
plot(f,20*log(X),f,20*log(orig_X),'--r');
box off;
grid on;
set(h,'Linewidth',1);
xlabel('Frequency(Hz)');
ylabel('Amplitude(dB)');
legend('3 Averages','Original');
title(['Frequency Spectrum w/ Noise Averaging: 3 Averages']);
axis([0 250 -50 120]);

%% Example of Noise Averaging, 10 Coherent Averages
fs = 500; % sample freq
Ts=1/fs; % sample period
n=0:1*fs; % time window (fs=samples per second), 1 second
tn=n*Ts; % sample index
x=cos(2*pi*50*tn);
k=.2; %Noise amplitude
n1=randn(size(x));
n2=randn(size(x));
n3=randn(size(x));
n4=randn(size(x));
n5=randn(size(x));
n6=randn(size(x));
n7=randn(size(x));
n8=randn(size(x));
n9=randn(size(x));
n10=randn(size(x));
s1=x+(k*n1);
s2=x+(k*n2);
s3=x+(k*n3);
s4=x+(k*n1);
s5=x+(k*n2);
s6=x+(k*n3);
s7=x+(k*n1);
s8=x+(k*n2);
s9=x+(k*n3);
s10=x+(k*n1);
average_s=((s1+s2+s3+s4+s5+s6+s7+s8+s9+s10)/10);
% Time domain plot of x(t)=cos(2*pi*50*tn)
figure('Color',[1 1 1]);
h=plot(tn,average_s,'g','LineWidth',2);box off;hold on;plot(tn,x,'--r');
title('Sinusoid w/ Additive White Noise: 10 Averages');
xlabel('Time(s)');
ylabel('Amplitude');
legend('composite signal','averaged');
% Frequency Domain Plot
figure('Color',[1 1 1]);
[X,f]=ComputeSpectrum(average_s,fs,2^12);
plot(f,20*log(X),f,20*log(orig_X),'--r');
box off;
grid on;
set(h,'Color',[0 0.8 1]);
set(h,'Linewidth',1);
xlabel('Frequency(Hz)');
ylabel('Amplitude(dB)');
title(['Frequency Spectrum w/ Noise Averaging: 10 Averages']);
legend('10 Averages','Original');
axis([0 250 -50 120]);

DSP Tutorial #5: Frequency Analysis

This example computes the spectrum of a complex signal using the ComputeSpectrum function.

Couple important things to notice:
- Our range of frequencies on the x-axis is Fs/2
- As Fs increases, the resolution in the frequency domain is much better
- As our FFT bins (2^N) increases, the more bins we have (we can still resolve the freq content provided we meet the sample theorem, we have more bins - basically a finer plot)
- As our windows length increases, the better resolution we have


close all;
clear all;
fs=500;
Ts=1/fs;
n=0:1*fs;
tn=n*Ts;
x=5+cos(2*pi*200*tn)+10+cos(2*pi*100*tn);
% Time Domain Plot
figure('Color',[1 1 1]);
plot(tn,x);box off;
title(['Time Domain Signal']);
[X,f]=ComputeSpectrum(x,fs,2^14);
figure('Color',[1 1 1]);
plot(f,20*log(X));
title(['Freq Domain']);
axis('tight');

The Compute Spectrum Function

This function uses the FFT algorithm to compute the spectrum of an input with a number of points:

function [X,f]=ComputeSpectrum(x,fs,N);
% ComputeSpectrum Compute the Spectrum
%
% x input signal
% fs Sampling frequency of the input signal
% N number of points to evaluate spectrum (integer number of 2)
%
% X Absolute value of the spectrum (capital X for spectrum)
% f Frequency axis of the spectrum
%
% This function uses the FFT algorithm to compute the spectrum of an
% input signal with a number of points (N) N-point FFT where N must be
% a power of 2 (e.g. 2^10)
%
%Example: Compute the spectrum of sinusoid at 3Hz, fs=36Hz
%
% x=cos(2*pi*3*[0:12*3]*1/(12*3)); 12 times max freq
% [X,f]=ComputeSpectrum(x,12*3,2^12);
%
%Revision:0.0.1 jeastham
%
X=abs(fft(x,N));
f=(fs/N)*[0:N-1]; %make the axis freq
X=X(1:end/2);
f=f(1:end/2);

DSP Tutorial #4: Plotting Complex Signals

In this example we will plot a more complex time domain signal in Matlab:



clear all; %clear all variables
fs = 50; % sample freq
Ts=1/fs; % sample period
n=0:1*fs; % time window (fs=samples per second), 1 second
tn=n*Ts; % sample index
sig_a=cos(2*pi*5*tn);
sig_b=1.5+cos(2*pi*20*tn);
x=sig_a+sig_b;
% Time domain plot of x(t)=cos(2*pi*5*tn)
figure('Color',[1 1 1]);
h=plot(tn,x,'g','LineWidth',3);box off;hold on;
plot(tn,sig_a,'--r','LineWidth',2);
plot(tn,sig_b,'--k','LineWidth',2);
title('Example of Composite Signal: 2 Signals Added');
legend('Composite','Signal 1','Signal 2');
xlabel('Time(s)');
ylabel('Amplitude');
figure('Color',[1 1 1]);
plot(tn,sig_a);
title('Stem Plot Example: Fs = 50 (50 Points)');
axis('tight');

Sum of Sinusoids Triangle Approximation

Triangle Approximation: Sum of Sinusoids Example

The following Matlab code provides an example of the fact that any signal can constructed using a sum of sinusoids.



clear all; %clear all variables
fs = 50; % sample freq
Ts=1/fs; % sample period
n=0:1*fs; % time window (fs=samples per second), 1 second
tn=n*Ts; % sample index
wave_1=sin(2*pi*5*tn);
wave_2=sin(-1/9*sin(6*pi*5*tn));
wave_3=sin(1/25*sin(10*pi*5*tn));
wave_4=sin(-1/49*sin(14*pi*5*tn));
wave_5=sin(1/81*sin(18*pi*5*tn));
wave_6=sin(-1/121*sin(22*pi*5*tn));
wave_7=sin(1/169*sin(26*pi*5*tn));
x = (8/(pi^2))*(wave_1+wave_2+wave_3+wave_4+wave_5+wave_6+wave_7);
% Time domain plot of x(t)=cos(2*pi*5*tn)
figure('Color',[1 1 1]);
h=plot(tn,wave_1,'--r',tn,wave_2,'--k',tn,x,'g','LineWidth',2);box off;hold on;
xlabel('Time(s)');
ylabel('Amplitude');
title('Square Wave Approximation');
axis([0 1 -1.5 1.5]);

Plotting a Discrete Time-Domain Signal in Matlab

In order to properly construct a sinusoid in Matlab we need to do a few things. We need to give Matlab the equation of our sinusoid. We also need to tell Matlab how often to sample the sinusoid and for how long we want to sample our sinusoid. Essentially, we want to tell Matlab “ok, watch this 5Hz sinusoid for 1 second. In that one second, take 50 discrete samples (actually 51), and plot the data. I should see a point every 0.02 seconds, for a total of 51 total points (Matlab starts at zero).” What we want to do may seem pretty simple but this process often confuses students. This overview will help in understanding the proper construction of sinusoids. I say “proper” because we will use the proper DSP terms when we tell Matlab what to do. For example, we will use terms like “sample frequency” and “sample period” to describe our signal.

To establish our signal we will need to:

1. Set the sample frequency FS
2. Set the sample period Ts
3. Set the time window (number of samples)
4. Set the sample index Tn (time index)
5. Establish the signal x(t) and plot
   
   

1. Set the Sample Frequency (Fs):

This is our sample frequency. This is basically how many “points” we will use to construct each signal. Remember it’s zero based, so it’s really plus one. In the example below, Fs was set to 50, so we have 51 points in our plot, we ‘sample’ the signal 51 times. This is a 5Hz signal, so we are sampling it 10 times per cycle, which satisfies Nyquist-Shannon’s sampling theorem.

Matlab code:

Fs = 50;

2. Set the Sample Period (Ts):

Next, we will set the sample period, Ts. The sample period is calculated as 1/Fs, or 1/50 in this example = 0.02. This is essentially the ‘granularity’ of our sampling. In other words we will have a discrete point every 0.02 units (which will be time). If we zoom into our signal and only look at the first 0.2 seconds, we can see that we have 11 datapoints, one data point every 0.02s (remember it’s zero based).

Matlab code:
Ts = 1/Fs;

3. Set the Time Window (n):

Next we will establish out how many points are in one second of our signal. We will setup a matrix from zero to one second. In our example Fs = 50, so n = 0:1*Fs. When we set n in Matlab, all we are really doing is setting up a matrix of numbers from 0 to 50. That is: 0,1,2,3,4,5,6,…..,50. We will use this matrix to set the sample index Tn in the next step.

n is now equal to the matrix below:

0 1 2 3 4 5 6 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 46 47 48 49 50

Matlab code:
n =0:1*Fs;

4. Establish the Sample Index (Tn):

Now we will tell Matlab how to chop up or signal. We set the Sample Index as Tn = Ts * n. In our example, Tn will equal .02. We want to see a discrete point every 0.02 seconds. We will take n, or matrix from 0 to 500 and multiply it by Ts which was 0.02. Now we have a matrix with 51 elements from 0 to 1 spaced every 0.02 seconds. We are now ready to plot our sinusoid.

Tn is now equal to the matrix below:

0 0.0200 0.0400 0.0600 0.0800 0.1000 0.1200 0.1400 0.1600 0.1800 0.2000 0.2200 0.2400 0.2600 0.2800 0.3000 0.3200 0.3400 0.3600 0.3800 0.4000 0.4200 0.4400 0.4600 0.4800 0.5000 0.5200 0.5400 0.5600 0.5800 0.6000 0.6200 0.6400 0.6600 0.6800 0.7000 0.7200 0.7400 0.7600 0.7800 0.8000 0.8200 0.8400 0.8600 0.8800 0.9000 0.9200 0.9400 0.9600 0.9800 1.0000

Matlab code:

Tn =Ts*n;

5. Plot our 5Hz Sinusoid:

Plotting our sinusoid is pretty easy at this point. Remember, we want a 5Hz cosine signal constructed of 51 points in one second; we want to “watch” the signal for 1 second and plot the results. We can plot the X axis correctly with time also using tn.

clear all; %clear all variables

fs = 50; % sample freq

Ts=1/fs; % sample period

n=0:1*fs; % time window (fs=samples per second), 1 second

tn=n*Ts; % sample index

sig_a=cos(2*pi*5*tn);
figure('Color',[1 1 1]);

stem(tn,sig_a);

title('Stem Plot Example: Fs = 50 (50 Points)');

axis('tight');

DSP Tutorial Matlab Script

The script below will get you started plotting basic time domain signals. The following tutorials are available on additional DSP topics:

• DSP Tutorial #1: Basic Plots
• DSP Tutorial #2: Signal Creation
• DSP Tutorial #3: Sum of Sinusoids
• DSP Tutorial #4: Complex Signals
• DSP Tutorial #5: Frequency Analysis
• DSP Tutorial #6: Looping
• DSP Tutorial #7: Audio Processing
• DSP Tutorial #8: Spectrum Function
• DSP Tutorial #9: Basic digital filter
• DSP Tutorial #10: Windowing

This is a simple Matlab script which describes the basic plotting of time-domain signals.

%% DSP Tutorial Script
% James Eastham
% Member, IEEE
% Created on: 01/08/2010
% Revision: R1

%% Plotting a Basic Time Domain Signals
% Sinewave with a peak amplitude of 1 and a fequency of freq
clear all;
close all;
freq = 5; %5Hz
n=0:.01:.2; %number of discrete values on X axis
x = cos (2*pi*freq*n); %our cosine signal
figure('Color',[1 1 1]);
plot(n,x);
title('Example of Tim-Domain Signal');
xlabel('Time(s)');
ylabel('Amplitude');

%% Plotting a Basic Time Domain Signal Continued
figure('Color',[1 1 1]);
stem(n,x,'--rs','LineWidth',2,'MarkerEdgeColor','k','MarkerFaceColor','g','Markersize',10);
hold on;
plot(n,x);
title('Example of Time-Domain Signal Plot');
xlabel('Time(s)');
ylabel('Amplitude');


%% Plotting with the Subplot Feature
% subplot(m,n,p) breaks the figure window into an
% m-by-n matrix
figure('Color',[1 1 1]);
subplot(2,1,1);
plot(n,x);
title('Example of the Time-Domain Signal');
xlabel('Time(s)');
ylabel('Amplitude');
ns=1:length(n);
subplot(2,1,2);
stem(ns,x,'.r');
title('Example of the Descrete Points');
xlabel('Descrete Sample');
ylabel('Amplitude');
axis([1 21 -1.2 1.2]);
grid on;