1、设计出一套完整的系统，对信号进行频谱分析和滤波处理； 1产生一个连续信号，包含低频，中频，高频分量，对其进行采样，进行频谱分析，分别设计三种高通，低通，带通滤波器对信号进行滤波处理，观察滤波后信号的频谱。 2.采集一段含有噪音的语音信号（可以录制含有噪音的信号，或者录制语音后再加进噪音信号） ，对其进行采样和频谱分析，根据分析结果设计出一合适的滤波器滤除噪音信号。%写上标题 %设计低通滤波器： N,Wc=buttord() %估算得到 Butterworth 低通滤波器的最小阶数 N 和 3dB 截止频率 Wc a,b=butter(N,Wc); %设计 Butterworth 低通滤波器 h
2、,f=freqz(); %求数字低通滤波器的频率响应 figure(2); % 打开窗口 2 subplot(221); %图形显示分割窗口 plot(f,abs(h); %绘制 Butterworth 低通滤波器的幅频响应图 title(巴氏低通滤波器); grid; %绘制带网格的图像 sf=filter(a,b,s); %叠加函数 S 经过低通滤波器以后的新函数 subplot(222); plot(t,sf); %绘制叠加函数 S 经过低通滤波器以后的时域图形 xlabel(时间 (seconds); ylabel(时间按幅度 ); SF=fft(sf,256); %对叠加函数 S 经
3、过低通滤波器以后的新函数进行 256 点的基2 快速傅立叶变换 w= %新信号角频率 subplot(223); plot(); %绘制叠加函数 S 经过低通滤波器以后的频谱图 title(低通滤波后的频谱图); %设计高通滤波器 N,Wc=buttord() %估算得到 Butterworth 高通滤波器的最小阶数 N 和 3dB 截止频率 Wc a,b=butter(N,Wc,high); %设计 Butterworth 高通滤波器 h,f=freqz(); %求数字高通滤波器的频率响应 figure(3); subplot(221); plot(); %绘制 Butterworth 高通
4、滤波器的幅频响应图 title(巴氏高通滤波器); grid; %绘制带网格的图像 sf=filter(); %叠加函数 S 经过高通滤波器以后的新函数 subplot(222); plot(t,sf); ;%绘制叠加函数 S 经过高通滤波器以后的时域图形 xlabel(Time(seconds); ylabel(Time waveform); w; %新信号角频率 subplot(223); plot(); %绘制叠加函数 S 经过高通滤波器以后的频谱图 title(高通滤波后的频谱图); %设计带通滤波器 N,Wc=buttord() %估算得到 Butterworth 带通滤波器的最小阶
5、数 N 和 3dB 截止频率 Wc a,b=butter(N,Wc); %设计 Butterworth 带通滤波器 h,f=freqz(); %求数字带通滤波器的频率响应 figure(4); subplot(221); plot(f,abs(h); %绘制 Butterworth 带通滤波器的幅频响应图 title(butter bandpass filter); grid; %绘制带网格的图像 sf=filter(a,b,s); %叠加函数 S 经过带通滤波器以后的新函数 subplot(222); plot(t,sf); %绘制叠加函数 S 经过带通滤波器以后的时域图形 xlabel(T
6、ime(seconds); ylabel(Time waveform); SF=fft(); %对叠加函数 S 经过带通滤波器以后的新函数进行 256 点的基2 快速傅立叶变换 w=( %新信号角频率 subplot(223); plot(); %绘制叠加函数 S 经过带通滤波器以后的频谱图 title(带通滤波后的频谱图);matlab 如何做频谱分析% -% 方案 1：“x = a*cos(2*pi*w*t)”的形式：% -% 注意：1.时域的持续时间范围应较大；% 2.频率 w 与序列 k 的对应关系（N 为序列总长度）：w = 1/dt * k/N;% 3.采样频率 1/dt 应大于
7、w 的 2 倍% 4.结果曲线的峰值的横坐标对应的就是 w 和-w 值% -a = 0.75;w = pi/3;dt = 0.2;t = -30*pi:dt:30*pi;N = size(t, 2);x = a*cos(2*pi*w*t);y = fft(x);y = fftshift(y);figure;subplot(2,1,1);plot(t,x);subplot(2,1,2);plot(1/dt*(-N/2+1:N/2)/N, abs(y);% -% 方案 2：“x = a*cos(w*t)”的形式：% -% 注意：1.时域的持续时间范围应较大；% 2.频率 w 与序列 k 的对应关系
8、（N 为序列总长度）：w = 1/dt *2*pi* k/N;% 3.采样频率 1/dt 应大于 w/(2*pi) 的 2 倍% 4.结果曲线的峰值的横坐标对应的就是 w 和-w 值% -a = 0.75;w = pi/3;dt = 1;t = -20*pi:dt:20*pi;N = size(t, 2);x = a*cos(w*t);y = fft(x);y = fftshift(y);figure;subplot(2,1,1);plot(t,x);subplot(2,1,2);plot(1/dt*2*pi*(-N/2:N/2-1)/N, abs(y);备注：由于使用了 fftshift，所
9、以得到的 频谱序列关于原点对称，如果不需要负半轴的话自行修改一下就 ok 了function f=frequency(x,fs)dtlen=length(x);t=(0:dtlen-1)/fs;subplot(211);plot(t,x);axis tight;y=abs(fft(x)*2/dtlen;ff=(0:dtlen/2-1)*fs/dtlen;subplot(212);f=y(1:floor(dtlen/2);plot(ff,f);axis tight;return;Matlab 编程实现 FFT 实践及频谱分析内容1用 Matlab 产生正弦波,矩形波,以及白噪声信号，并显示各自时
10、域波形图2进行 FFT 变换，显示各自频谱图，其中采样率，频率、数据长度自选3做出上述三种信号的均方根图谱,功率图谱,以及对数均方根图谱4用 IFFT 傅立叶反变换恢复信号，并显示恢复的正弦信号时域波形图源程序%*% FFT 实践及频谱分析 %*%*%*1.正弦波*%fs=100;%设定采样频率N=128;n=0:N-1;t=n/fs;f0=10;%设定正弦信号频率%生成正弦信号x=sin(2*pi*f0*t);figure(1);subplot(231);plot(t,x);%作正弦信号的时域波形xlabel(t);ylabel(y);title(正弦信号 y=2*pi*10t 时域波形);
11、grid;%进行 FFT 变换并做频谱图y=fft(x,N);%进行 fft 变换mag=abs(y);%求幅值f=(0:length(y)-1)*fs/length(y);%进行对应的频率转换figure(1);subplot(232);plot(f,mag);%做频谱图axis(0,100,0,80);xlabel(频率(Hz);ylabel(幅值);title(正弦信号 y=2*pi*10t 幅频谱图 N=128);grid;%求均方根谱sq=abs(y);figure(1);subplot(233);plot(f,sq);xlabel(频率(Hz);ylabel(均方根谱);title
12、(正弦信号 y=2*pi*10t 均方根谱);grid;%求功率谱power=sq.2;figure(1);subplot(234);plot(f,power);xlabel(频率(Hz);ylabel(功率谱);title(正弦信号 y=2*pi*10t 功率谱);grid;%求对数谱ln=log(sq);figure(1);subplot(235);plot(f,ln);xlabel(频率(Hz);ylabel(对数谱);title(正弦信号 y=2*pi*10t 对数谱);grid;%用 IFFT 恢复原始信号xifft=ifft(y);magx=real(xifft);ti=0:len
13、gth(xifft)-1/fs;figure(1);subplot(236);plot(ti,magx);xlabel(t);ylabel(y);title(通过 IFFT 转换的正弦信号波形);grid;%*2.矩形波*%fs=10;%设定采样频率t=-5:0.1:5;x=rectpuls(t,2);x=x(1:99);figure(2);subplot(231);plot(t(1:99),x);%作矩形波的时域波形xlabel(t);ylabel(y);title(矩形波时域波形);grid;%进行 FFT 变换并做频谱图y=fft(x);%进行 fft 变换mag=abs(y);%求幅值
14、f=(0:length(y)-1)*fs/length(y);%进行对应的频率转换figure(2);subplot(232);plot(f,mag);%做频谱图xlabel(频率(Hz);ylabel(幅值);title(矩形波幅频谱图);grid;%求均方根谱sq=abs(y);figure(2);subplot(233);plot(f,sq);xlabel(频率(Hz);ylabel(均方根谱);title(矩形波均方根谱);grid;%求功率谱power=sq.2;figure(2);subplot(234);plot(f,power);xlabel(频率(Hz);ylabel(功率谱
15、);title(矩形波功率谱);grid;%求对数谱ln=log(sq);figure(2);subplot(235);plot(f,ln);xlabel(频率(Hz);ylabel(对数谱);title(矩形波对数谱);grid;%用 IFFT 恢复原始信号xifft=ifft(y);magx=real(xifft);ti=0:length(xifft)-1/fs;figure(2);subplot(236);plot(ti,magx);xlabel(t);ylabel(y);title(通过 IFFT 转换的矩形波波形);grid;%*3.白噪声*%fs=10;%设定采样频率t=-5:0.
16、1:5;x=zeros(1,100);x(50)=100000;figure(3);subplot(231);plot(t(1:100),x);%作白噪声的时域波形xlabel(t);ylabel(y);title(白噪声时域波形);grid;%进行 FFT 变换并做频谱图y=fft(x);%进行 fft 变换mag=abs(y);%求幅值f=(0:length(y)-1)*fs/length(y);%进行对应的频率转换figure(3);subplot(232);plot(f,mag);%做频谱图xlabel(频率(Hz);ylabel(幅值);title(白噪声幅频谱图);grid;%求均
17、方根谱sq=abs(y);figure(3);subplot(233);plot(f,sq);xlabel(频率(Hz);ylabel(均方根谱);title(白噪声均方根谱);grid;%求功率谱power=sq.2;figure(3);subplot(234);plot(f,power);xlabel(频率(Hz);ylabel(功率谱);title(白噪声功率谱);grid;%求对数谱ln=log(sq);figure(3);subplot(235);plot(f,ln);xlabel(频率(Hz);ylabel(对数谱);title(白噪声对数谱);grid;%用 IFFT 恢复原始信号xifft=ifft(y);magx=real(xifft);ti=0:length(xifft)-1/fs;figure(3);subplot(236);plot(ti,magx);xlabel(t);ylabel(y);title(通过 IFFT 转换的白噪声波形);grid;
