# Post

## Done :)

That should be everything you need to extract a sinusoid!

## MATLAB Simulation Code

The following MATLAB code should generate and apply your filter:

clear all; close all; clc;

s = tf('s');       %declare s as Laplace variable

%all times in seconds - freq in Hz.
tmax = 10;         %time window of simulation
dt = 0.05;         %dt of simulation
fsquare = 1;       %frequency of square wave
Amp = 1;           %amplitude (peak to peak)

wc = 2*pi*1;       %cutoff frequency of filter
zeta = 0.2 ;       %damping const. of filter

fftsize = 2^16;    %number of points the generates

G_filt = 1/( 1   +   2*zeta*s/wc   +   s^2/(wc^2) );

%make a square wave
t = linspace(0,tmax,tmax/dt);
x = 1:(tmax/dt);
k = 1/(2*dt*fsquare);
y_square = (Amp*(mod(x,2*k)- mod(x,k))/k - 0.5*Amp);

%get filter's impulse repsonse
y_impulse = impulse(G_filt,t);

%apply filter (convolve then truncate)
y_filtered = conv(y_square,y_impulse)*dt;
y_filtered = y_filtered(1:length(y_square));

%plot bode response of filter
figure();
bode(G_filt);
grid; grid minor;

%plot impulse response of filter
figure();
plot(t,y_impulse,'g','linewidth',2);
title('Filter impulse response')
xlabel('Time')
ylabel('Intensity / Volts')
grid; grid minor;

%plot original signal
figure();
plot(t,y_square,'m','linewidth',2);
xlabel('Time')
ylabel('Intensity / Volts')
hold on

%plot filtered response
plot(t,y_filtered,'b','linewidth',2);
xlabel('Time')
ylabel('Intensity / Volts')
grid; grid minor;
legend('Original Square','Filtered Signal','Location','southwest')

%make user pick point at which transients are gone
uiwait(msgbox('Pick the point on the graph where the transients are gone.'));
[t_truncate,~] = ginput(1);
n_truncate = round(t_truncate/dt);

%truncate the filtered dataset and shift to zero to remove DC component.
y_filtered_truncated = y_filtered(n_truncate:end);
t_truncated = t(n_truncate:end);
y_dc = trapz(y_filtered_truncated,t_truncated)/(t_truncated(end) - t_truncated(1));

%run ffts
Y_fil = fft(y_filtered_truncated,fftsize);
Y_sq = fft(y_square(n_truncate:end),fftsize);

%evaluate power spectra
Pyy_sq = Y_sq.*conj(Y_sq)/fftsize;
Pyy_fil = Y_fil.*conj(Y_fil)/fftsize;
freq = (1/dt)*(0:(fftsize/2))/fftsize;

%plot fft for square
figure();
plot(freq,Pyy_sq(1:(1+fftsize/2)),'k','LineWidth', 2);
xlabel('Frequency (Hz)');
ylabel('Power');
hold on

%plot fft for filtered data
plot(freq,Pyy_fil(1:(1+fftsize/2)),'r','LineWidth', 2);
title('Signal FFTs')
xlabel('Frequency (Hz)');
ylabel('Power');
grid; grid minor;
legend('Square Wave','Filtered Signal','Location','northeast')
* This program runs in MATLAB 2017b

1. MATLAB code available at bottom of this page ^ right here! ↩︎