http://s13/mw690/0035H5bGzy6ZDN0IIiMfc&690
http://s5/mw690/0035H5bGzy6ZDN3Gq6U14&690

%频率不同的三个正弦波信号相加，分别输出时域波形和频域频谱 three_cos.m

%设计低通滤波器，改变低通滤波器截止频率，输出不同波形

%1.波形产生程序 three_cos(t)；

%2.傅里叶正变换[f, sf]=FFT_SHIFT(t, st)；

%3.傅里叶反变换[t,st]=IFFT_SHIFT(f,Sf)；

%4.低通滤波器[t,st]=RECT_LPF(f,Sf,B)；

 

% Signal

dt=0.01;    %时间采样间隔

T=5;        %信号持续时间

N=floor(T/dt);  %点数

t=[0:N-1]*dt;   %t

 

f1=1;f2=5;f3=10; %信号不同频率值

m1=cos(2*pi*f1*t);

m2=cos(2*pi*f2*t);

m3=cos(2*pi*f3*t);

mt1=m1+m2+m3;

[f,mt1f]=FFT_SHIFT(t,mt1); %2.傅里叶正变换 

%输出:f-抽样频率；mt1f-频率；输入:t-时间；mt1-输入波形

fmax=2;    

f=1;f2=5;f3=10; %信号不同频率值

B1=fmax %设置低通滤波器带宽，准备调用低通滤波器

[t,mt1_t]=RECT_LPF(f,mt1f,B1); %低通滤波器滤除高频，由B1决定

[f,mt1_tf]=FFT_SHIFT(t,mt1_t); %2.傅里叶正变换 

 

 %时域、频域显示

subplot(511);

plot(t,m1);

ylabel('m1(t)');

title('三个不同频率正弦波合成波和频谱');

subplot(512);

plot(t,m2);

ylabel('m2(t)');

subplot(513);

plot(t,m3);

ylabel('m3(t)');

subplot(514);

plot(t,mt1);

ylabel('m(t)');

subplot(515);

plot(f,mt1f);

ylabel('m(f)');

axis([-40,40,0,3]);

 

figure(2)

subplot(511);

plot(t,m1);

ylabel('m1(t)');

title('低通滤波器选频合成波和频谱');

subplot(512);

plot(t,m2);

ylabel('m2(t)');

subplot(513);

plot(t,m3);

ylabel('m3(t)');

subplot(514);

plot(t,mt1_t);

ylabel('m(t)');

subplot(515);

plot(f,mt1_tf);

ylabel('fmax=2');

axis([-40,40,0,3]);

%=======================================

% 调用函数  傅里叶变换

function [f, sf]=FFT_SHIFT(t, st)

%This function is FFT to calculate a signal’s Fourier transform

%Input: t: sampling time , st : signal data. Time length must greater than 2

%output: f : sampling frequency , sf: frequency

%output is the frequency and the signal spectrum

dt=t(2)-t(1);

T=t(end);

df=1/T;

N=length(t);

f=[-N/2:N/2-1]*df;

sf=fft(st);

sf=T/N*fftshift(sf);

 

% 傅里叶反变换

function  [t,st]=IFFT_SHIFT(f,Sf)

df=f(2)-f(1);

fmax=(f(end)-f(1)+df);

dt=1/fmax;

N=length(f);

t=[0:N-1]*dt;

Sf=fftshift(Sf);

st=fmax*ifft(Sf);

st=real(st);

 

% 赵鸿图-低通滤波器

function [t,st]=RECT_LPF(f,Sf,B)

df=f(2)-f(1);

fN=length(f);

RectH=zeros(1,fN);

BN=floor(B/df);

BN_SHIFT=[-BN:BN-1]+floor(fN/2);

RectH(BN_SHIFT)=1;

Yf=RectH.*Sf;

[t,st]=IFFT_SHIFT(f,Yf);

 

%参考书： 赵鸿图 茅艳   通信原理MATLAB仿真教材[M]. 人民邮电出版社