演習3解答(12 月 1 日まで) Q3.m --------------------------------------------------

演習3解答(12 月 1 日まで)
Q3.m -----------------------------------------------------------------------------% 孤立三角パルスのフーリエ変換
clear all
close all
df=0.02;
fs=10;
ts=1/fs;
t=-5:ts:5;
%
%
%
%
周波数分解能 [Hz]
標本化周波数 [Hz]
標本化間隔 [sec]
時間 [sec]
% 三角パルス x1(t)
x1=zeros(size(t));
x1(41:51)=(0:1:10)/10;
x1(51:60)=(10:-1:1)/10;
% 三角パルス x2(t)
x2=zeros(size(t));
x2(51:61)=(0:10)/10;
x2(61:70)=(10:-1:1)/10;
figure()
plot(t,x1,'-',t,x2,'--')
axis([-5 5 -1 2])
legend('x1(t)','x2(t)')
xlabel('sec')
% 三角パルス x1(t),x2(t)をプロット
% x軸とy軸の範囲を指定
% 凡例を指定
[X1,x11,df1]=fourier_analysis(x1,ts,df);
[X2,x21,df2]=fourier_analysis(x2,ts,df);
X11=fftshift(X1)/fs;
X21=fftshift(X2)/fs;
f=[0:df1:df1*(length(x11)-1)]-fs/2;
%
%
%
%
%
三角パルス x1(t)をMATLAB DFT(FFT)
三角パルス x2(t)をMATLAB DFT(FFT)
MATLAB DFT(FFT)とフーリエ変換の関係
MATLAB DFT(FFT)とフーリエ変換の関係
スペクトルに対応した周波数軸を作成
figure()
plot(f,abs(X11),'k-',f,abs(X21),'k--')
% 振幅スペクトル(実線と破線が重なる)
legend('x1(t)の振幅スペクトル','x2(t)の振幅スペクトル')
xlabel('Hz')
figure() % 位相スペクトル [位相シフトの確認]
plot(f(250:264),angle(X11(250:264)),'k-',f(250:264),angle(X21(250:264)),'k--')
legend('x1(t)の位相スペクトル','x2(t)の位相スペクトル')
axis([-0.15 0.15 -4 4])
xlabel('Hz')
fourier_analysis.m -----------------------------------------------------------------function [M,m,df]=fourier_analysis(m,ts,df)
% m: 入力信号/出力信号
% ts: 標本時間間隔
% df: 必要周波数分解能/実際の周波数分解能
% M: 周波数スペクトル
fs=1/ts;
% 入力信号の標本周波数
n1=fs/df;
% 必要周波数分解能を得るためのサンプル数
n2=length(m);
% 入力信号の長さ
n=2^max(nextpow2(n1),nextpow2(n2)); % 必要周波数分解能を得るためのfftサンプル数
M=fft(m,n);
% 周波数スペクトル
m=[m,zeros(1,n-n2)];
% 0をつけたデータ
df=fs/n;
% 実際の周波数分解能