Matlab link : https://drive.google.com/open?id=0B5MDnBF8A-aENF9YLTRZaXZxd3c The
ID: 3868465 • Letter: M
Question
Matlab link : https://drive.google.com/open?id=0B5MDnBF8A-aENF9YLTRZaXZxd3c
The first tab of the Excel data file for this problem contains a complex periodic function composed of three sine waves with different amplitudes and frequencies, sampled at 50 Hz for 20 seconds. Using the supplied MATLAB function MAE2381_Su17_HW4_FilterProgram.m import the data into MATLAB and then use a filter to eliminate the amplitude content at the lowest and highest frequencies. Only the amplitude content at the middle frequency should be left over once the filter is applied. Figures 1 and 2 show the complex periodic function and its frequency spectrum with a line that indicates the part of the signal that should be left over after filtering. a) Using the frequency spectrum plot of the original function (Fig. 2, or the same from MATLAB), write the equation of the of the complex periodic function and identify the frequencies of the three component sinusoids b) Using MAE2381_Sul7_HW4_FilterProgram.m, apply a filter to the data to isolate the middle sinusoid. The program is already completed and will generate the required figures automatically when run: you can read the comments to learn how it works. The filter is a k-order, Butterworth, bandpass type. Lines 54-56 contain filter settings that you must change in order to filter out the low- and high-frequency sinusoids. Once you believe you have correctly filtered the data, list what filter settings you used.Explanation / Answer
Fs = 150; % Sampling frequency
t = 0:1/Fs:1; % Time vector of 1 second
f = 5; % Create a sine wave of f Hz.
i (i*f)
-0.5
0
Amplitude
x = sin(2* pi*t*f);
nfft = 1024; % Length of FFT
% Take fft, padding with zeros so that length(X)
is equal to nfft
X = fft(x,nfft);
% second 0 0.2 0.4 0.6 0.8 1
-1
Time (s)
80
Power Spectrum of a Sine Wave
FFT is symmetric, throw away half
X = X(1:nfft/2);
% Take the magnitude of fft of x
mx = abs(X);
% Frequency vector
f = (0:nfft/2-1)*Fs/nfft;
40
60
Power
% Generate the plot, title and labels.
figure(1);
plot(t,x);
title('Sine Wave Signal');
xlabel('Time (s)');
l b l('A lit d ')
0 10 20 30 40 50 60 70 80
0
20
Frequency P
ylabel('Amplitude');
figure(2);
plot(f,mx);
title('Power Spectrum of a Sine Wave');
xlabel('Frequency (Hz)');
(Hz)
ylabel('Power');
Fs = 60; % Sampling frequency
t = -.5:1/Fs:.5;
x = 1/(sqrt(2*pi*0.01))*(exp(-t.^2/(2*0.01)));
nfft = % FFT
1
2
Amplitude
1024; Length of % Take fft, padding with zeros so that
length(X) is equal to nfft
X = fft(x,nfft);
% FFT is symmetric, throw away second half
X = X(1:nfft/2);
-0.5 0 0.5
0
Time (s)
( )
% Take the magnitude of fft of x
mx = abs(X);
% This is an evenly spaced frequency vector
f = (0:nfft/2-1)*Fs/nfft;
% Generate the plot, title and labels.
40
50
60
Power Spectrum of a Gaussian Pulse
r
figure(1);
plot(t,x);
title('Gaussian Pulse Signal');
xlabel('Time (s)');
ylabel('Amplitude');
figure(2);
0
10
20
30
Power
plot(f,mx);
title('Power Spectrum of a Gaussian Pulse');
xlabel('Frequency (Hz)');
ylabel('Power');
Fs = 200; % Sampling frequency
t = 0:1/Fs:1; % Time vector of 1 second
x = chirp(t,0,1,Fs/6);
nfft = 1024; % Length of FFT
-0.5
0
Amplitude
; g
% Take fft, padding with zeros so that
length(X) is equal to nfft
X = fft(x,nfft);
% FFT is symmetric, throw away second half
X = X(1:nfft/2);
0 0.2 0.4 0.6 0.8 1
-1
Time (s)
Power Spectrum of Chirp Signal
% Take the magnitude of fft of x
mx = abs(X);
% This is an evenly spaced frequency
vector
f = (0:nfft/2-1)*Fs/nfft;
% plot labels
15
20
25
p p g
wer
Generate the plot, title and labels.
figure(1);
plot(t,x);
title('Chirp Signal');
xlabel('Time (s)');
ylabel('Amplitude');
0 20 40 60 80 100
0
5
10
Pow
y p
figure(2);
plot(f,mx);
title('Power Spectrum of Chirp Signal');
xlabel('Frequency (Hz)');
ylabel('Power');
load quakevibration.mat
Fs = 1e3; % sample rate
NFFT = 512; % number of FFT points
segmentLength = 64; % segment length
% open loop acceleration power spectrum
[P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,'power');
% closed loop acceleration power spectrum
P1_CL = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,'power');
helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),...
'Frequency in Hz','Acceleration Power Spectrum in dB',...
'Resolution bandwidth = 15.625 Hz',{'Open loop', 'Closed loop'},[0 100])
NFFT = 512; % number of FFT points
segmentLength = 512; % segment length
[P1_OL,F] = pwelch(gfloor1OL,ones(segmentLength,1),0,NFFT,Fs,'power');
P1_CL = pwelch(gfloor1CL,ones(segmentLength,1),0,NFFT,Fs,'power');
helperFrequencyAnalysisPlot2(F,10*log10([(P1_OL) (P1_CL)]),...
'Frequency in Hz','Acceleration Power Spectrum in dB',...
'Resolution bandwidth = 1.95 Hz',{'Open loop', 'Closed loop'},[0 100])
% Take the magnitude of each FFT component of the signal
Yzp = abs(Y);
helperFrequencyAnalysisPlot1(F,abs(Yzp),unwrap(angle(Yzp)),NFFT,[],...
'Phase has been set to zero')
yzp = ifft(Yzp,'symmetric');
hplayer = audioplayer(yzp, Fs);
play(hplayer);
load officetemp.mat
Fs = 1/(60*30); % Sample rate is 1 sample every 30 minutes
t = (0:length(temp)-1)/Fs;
helperFrequencyAnalysisPlot2(t/(60*60*24*7),temp,...
'Time in weeks','Temperature (Fahrenheit)')
NFFT = length(temp); % Number of FFT points
F = (0 : 1/NFFT : 1/2-1/NFFT)*Fs; % Frequency vector
TEMP = fft(temp,NFFT);
TEMP(1) = 0; % remove the DC component for better visualization
helperFrequencyAnalysisPlot2(F*60*60*24*7,abs(TEMP(1:NFFT/2)),...
'Frequency (cycles/week)','Magnitude',[],[],[0 10])
load ampoutput1.mat
Fs = 3600;
NFFT = length(y);
% Power spectrum is computed when you pass a 'power' flag input
[P,F] = periodogram(y,[],NFFT,Fs,'power');
helperFrequencyAnalysisPlot2(F,10*log10(P),'Frequency in Hz',...
'Power spectrum (dBW)',[],[],[-0.5 200])
PdBW = 10*log10(P);
power_at_DC_dBW = PdBW(F==0) % dBW
[peakPowers_dBW, peakFreqIdx] = findpeaks(PdBW,'minpeakheight',-11);
peakFreqs_Hz = F(peakFreqIdx)
peakPowers_dBW
load ampoutput2.mat
SegmentLength = NFFT;
% Power spectrum is computed when you pass a 'power' flag input
[P,F] = pwelch(y,ones(SegmentLength,1),0,NFFT,Fs,'power');
helperFrequencyAnalysisPlot2(F,10*log10(P),'Frequency in Hz',...
'Power spectrum (dBW)',[],[],[-0.5 200])