firhamming_lpf.m

this a DSP program for FIR

%writen by Ziwei Chen 
clc;clear all;close all;
%产生模拟信号x(t)
t=[0:0.00001:1];
f1=200;f2=1800;f3=3000;
xt=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);
h1=figure;
set(h1,'name','进行滤波的输入信号');
subplot(2,1,1);
plot(t,xt);
xlabel('时间t （单位：s）');ylabel('幅度x(t)');
title('模拟时间信号x(t)');

%计算采样周期Ts及采样点数N,对原时间信号进行截断，选取Tp时间进行观察
fs=10000;Ts=1/fs;
%F=10;
N=128;%fs/F;                               %N=input('请输入采样点数N＝');   
Tp=N*Ts;                              %最小记录时间（对信号的观察时间）

%对模拟信号进行采样，得到序列信号x(n)
n=[0:N-1];
xn=sin(2*pi*f1*n*Ts)+sin(2*pi*f2*n*Ts)+sin(2*pi*f3*n*Ts);
subplot(2,1,2);
stem(n,xn);axis([0 1024 -3 3]);
xlabel('时间n');ylabel('幅度x(n)');
title('离散时间信号x(n)')

%对序列信号x(n)进行谱分析，得到信号的频谱
K=[0:N-1];
Nfft=N;
XK=fft(xn,Nfft);
h2=figure;
set(h2,'name','对输入信号进行频谱分析的结果（低通滤波前）');
plot(fs*K/Nfft,abs(XK));grid on;
set(gca,'XTickMode','manual','XTick',[100,1400,3000,4000,5000,6600,7900]);
xlabel('频率f  (单位：Hz)');ylabel('幅度|X(j2\pif)|');
title('输入信号的幅度谱（离散频谱）');
%----------------------------------------------------------------

%  Lowpass filter design - Hamming window
wp = 2*pi*250/fs;               %fpass=375Hz wp=0.075*pi
ws = 2*pi*1000/fs;               %fstop=1000Hz  deltaf=625Hz  ws=0.2*pi 

tr_width = ws - wp               %tr_width =0.125*pi

M = ceil(6.6*pi/tr_width) + 1

n=[0:1:M-1];

wc = (ws+wp)/2

hd = ideallp(wc,M);

w_ham = (hamming(M))';

h = hd .* w_ham;

[db,mag,pha,grd,w] = freqz_m(h,[1]);

delta_w = 2*pi/1000;

Rp = -(min(db(1:1:wp/delta_w+1)))        % Passband Ripple

As = -round(max(db(ws/delta_w+1:1:501))) % Min Stopband attenuation

% plots
h3=figure;
set(h3,'name','数字低通滤波器设计（采用哈明窗）');
subplot(1,1,1)

subplot(2,2,1); stem(n,hd); title('Ideal Impulse Response')

axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('hd(n)')

subplot(2,2,2); stem(n,w_ham);title('Hamming Window')

axis([0 M-1 0 1.1]); xlabel('n'); ylabel('w(n)')

subplot(2,2,3); stem(n,h);title('Actual Impulse Response')

axis([0 M-1 -0.1 0.3]); xlabel('n'); ylabel('h(n)')

subplot(2,2,4); plot(w/pi,db);title('Magnitude Response in dB');grid

axis([0 1 -100 10]); xlabel('frequency in pi units'); ylabel('Decibels')

set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1])

set(gca,'YTickMode','manual','YTick',[-50,0])

set(gca,'YTickLabelMode','manual','YTickLabels',['50';' 0'])
%----------------------------------------------------------------

%对输入序列进行低通滤波
HK=fft(h,Nfft);
YK=XK.*HK;
%y=xcorr(xn,ifft(H,Nfft));
%YK=fft(y,Nfft);
h4=figure;
set(h4,'name','对输入信号进行低通滤波的结果（低通滤波后）');
plot(fs*K/Nfft,abs(YK));grid on;%axis tight;
set(gca,'XTickMode','manual','XTick',[100,1400,3000,4000,5000,6600,7900]);
xlabel('频率f  (单位：Hz)');ylabel('幅度|Y(j2\pif)|');
title('输出信号的幅度谱（离散频谱）');