数字信号处理快速傅离叶变换和在信号处理中的应用实验

最新推荐文章于 2024-01-06 01:38:05 发布

小王yue

最新推荐文章于 2024-01-06 01:38:05 发布

阅读量2.1k

点赞数 1

分类专栏：音频处理 matlab 文章标签： matlab fft

本文链接：https://blog.csdn.net/qq_51940383/article/details/122512680

版权

音频处理同时被 2 个专栏收录

3 篇文章 0 订阅

订阅专栏

matlab

1 篇文章 0 订阅

订阅专栏

一、实验目的

学习和掌握基-2DIT-FFT算法原理和流图特点；
掌握根据算法特点实现FFT的方法；
对一段语音信号进行FFT变换，分别提取高频成分和低频成分，掌握采用FFT进行信号处理的初步方法。

二、实验内容

根据FFT流图特点编程实现FFT和IFFT（不允许使用MATLAB自带函数），画出幅值谱和相位谱；
比较直接计算DFT方法和FFT方法的结果，研究二种方法的运算速度
编程实现利用重叠相加法计算长序列和短序列的卷积。
对一段语音信号进行处理，分别提取高频和低频成分，观察提取出来信号的幅值谱和时域波形，分析产生的变化和听起来的感觉。

三、实验报告要求

1．步骤完整，按照报告已有项目填写，实验前需预习，掌握实验内容的原理。

2．实验报告必须包括：程序流程图、自己选择例子看程序运行结果进行分析。

（1）自己编制对某一个有限长度为N的序列直接计算其离散傅立叶变换的基-2的DIT-FFT程序，运行后画出幅值谱和相位谱图形；

（2）将上面程序编成子函数的形式，编制利用FFT来计算IFFT子程序；

（3）利用重叠相加法计算长序列和短序列的卷积，观察与不分段的用时差异。

（4）对一段语音信号进行处理，分别提取高频和低频成分，观察提取出来信号的幅值谱和时域波形，分析产生的变化和听起来的感觉，分析音乐信号本身的频率特点，变化不同截止频率详细分析所得结果。

（5）如果对语音信号不断下采样，幅值谱会发生什么样的变化，语音信号在时域和听觉感受上会发生什么的变化，分析出原因。

四、实验内容

图1 FFT和IFFT得到的幅值谱和相位谱

图2 DFT和FFT得到的幅值谱和相位谱

图3 DFT和FFT分别运行的时间

图4 不分段卷积与重叠相加法得到的结果

图5 不分段卷积与重叠相加法使用的时间

图6 初始语音信号时域频域图

图7 语音信号加噪后时域频域图

图8 语音信号低通滤波前后时域波形图

图9 语音信号低通滤波前后频谱图

图10 语音信号高通滤波前后时域波形图

图11 语音信号高通滤波前后频谱图

（1）fft代码

function y=ditfft(x,x_length)
%本程序对输入序列实现DIT-FFT基2算法，点数取大于等于长度的2的幂次
m=nextpow2(x_length);                       %求的x长度对应的2的最低幂次m[1,2,3,4,4,3,2,1,4,5],10
N=2^m;
if length(x)<N
    x=[x,zeros(1,N-length(x))];              %若的长度不是2的幂，补0到2的整数幂
end
nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1;   %求1：2^m数列的倒序
y=x(nxd);                                    %将倒序排列作为的初始值
for L=1:m                                    %将DFT做m次基2分解，从左到右，对每次分解作DFT运算
    D=2^L;
    u=1;                                     %旋转因子u初始化
    WN=exp(-1i*2*pi/D);                      %本次分解的基本DFT因子WN＝exp(-i*2*pi/Nmr)
      for j=1:D/2                            %本次跨越间隔内的各次碟形运算
          for k=j:D:N                        %本次碟形运算的跨越间隔为Nmr=2^mm
              kp=k+D/2;                      %确定碟形运算的对应单元下标
              t=y(kp)*u;                     %碟形运算的乘积项
              y(kp)=y(k)-t;                  %碟形运算的加法项
              y(k)=y(k)+t;
          end
          u=u*WN;                            %修改旋转因子，多乘一个基本DFT因子WN
      end
end
y=y(1:x_length);

（2）ifft代码

function y=ditifft(x,x_length)
%本程序对输入序列实现DIT-FFT基2算法，点数取大于等于长度的2的幂次
m=nextpow2(x_length);                       %求的x长度对应的2的最低幂次m
N=2^m;
if length(x)<N
    x=[x,zeros(1,N-length(x))];              %若的长度不是2的幂，补0到2的整数幂
end
nxd=bin2dec(fliplr(dec2bin([1:N]-1,m)))+1;   %求1：2^m数列的倒序
y=x(nxd);                                    %将倒序排列作为的初始值
for L=1:m                                    %将DFT做m次基2分解，从左到右，对每次分解作DFT运算
    D=2^L;
    u=1;                                     %旋转因子u初始化
    WN=exp(i*2*pi/D);                      %本次分解的基本DFT因子WN＝exp(i*2*pi/Nmr)
      for j=1:D/2                            %本次跨越间隔内的各次碟形运算
          for k=j:D:N                        %本次碟形运算的跨越间隔为Nmr=2^mm
              kp=k+D/2;                      %确定碟形运算的对应单元下标
              t=y(kp)*u;                     %碟形运算的乘积项
              y(kp)=y(k)-t;                  %碟形运算的加法项
              y(k)=y(k)+t;
          end
          u=u*WN;                            %修改旋转因子，多乘一个基本DFT因子WN
      end
end
y=y(1:x_length)/N;

（3）重叠相加法

function [Y] = overpl(x,h,N)
%N为分段段长度；x为长序列；h为短序列
Lx=length(x);
M=length(h);
x=[x,zeros(1,N)];
t=zeros(1,M-1);
Y=zeros(1,Lx+M-1);
a=floor(Lx/N);
y2=ditfft(h,N+M-1);
for K=0:a   %循环a+1次
    A=x(K*N+1:K*N+N);
    y1=ditfft(A,N+M-1);
    y3=y1.*y2;
    q=ditifft(y3,N+M-1);
    Y(K*N+1:K*N+M-1)=q(1:M-1)+t(1:M-1);
    Y(K*N+M:K*N+N)=q(M:N);
    t(1:M-1)=q(N+1:N+M-1);
end
Y=Y(1:Lx+M-1);

（4）测试代码

N=64;
n=0:N-1;
x1n=0.8*n;
% figure(1);
% subplot(2,3,1);stem(n,abs(x1n),'.');title("原始信号时域");xlabel('时间轴');ylabel('幅值')
% subplot(2,3,2);stem(n,angle(x1n),'.');title("原始信号频域");xlabel('时间轴');ylabel('频率幅值')
% X1k=ditfft(x1n,N);
% subplot(2,3,3);stem(n,abs(X1k),'.');title("ditfft时域");xlabel('时间轴');ylabel('幅值')
% subplot(2,3,4);stem(n,angle(X1k),'.');title("ditfft频域");xlabel('时间轴');ylabel('频率幅值')
% x1n1=ditifft(X1k,N);
% subplot(2,3,5);stem(n,abs(x1n1),'.');title("ditifft时域");xlabel('时间轴');ylabel('幅值')
% subplot(2,3,6);stem(n,angle(x1n1),'.');title("ditifft频域");xlabel('时间轴');ylabel('频率幅值')

N1=6;
n1=0:N1-1;
x2n=0.3*n1;
tic
Y=juanji(x1n,x2n);%卷积
toc
figure(2);
subplot(2,2,1);stem(0:N+N1-2,abs(Y),'.');title("卷积时域");xlabel('时间轴');ylabel('幅值')
subplot(2,2,2);stem(0:N+N1-2,angle(Y),'.');title("卷积频域");xlabel('时间轴');ylabel('频率幅值')
tic
Y1=overpl(x1n,x2n,N1);%重叠相加法
toc
subplot(2,2,3);stem(0:N+N1-2,abs(Y1),'.');title("重叠相加法时域");xlabel('时间轴');ylabel('幅值')
subplot(2,2,4);stem(0:N+N1-2,angle(Y1),'.');title("重叠相加法频域");xlabel('时间轴');ylabel('频率幅值')

figure(3);
tic
X2k=dft(x1n,N);
toc
subplot(2,2,1);stem(n,abs(X2k),'.');title("dft时域");xlabel('时间轴');ylabel('幅值')
subplot(2,2,2);stem(n,angle(X2k),'.');title("dft频域");xlabel('时间轴');ylabel('频率幅值')
tic
X1k=ditfft(x1n,N);
toc
subplot(2,2,3);stem(n,abs(X2k),'.');title("ditfft时域");xlabel('时间轴');ylabel('幅值')
subplot(2,2,4);stem(n,angle(X2k),'.');title("ditfft频域");xlabel('时间轴');ylabel('频率幅值')

（5）原始语音信号的时域频域波形代码

[x,fs]=audioread('D:\study\大三上\数字信号处理\实验\二\data1.wav');
n=length(x);
y=fft(x,n);
f=fs*(0:n/2-1)/n;
figure(1)subplot(2,1,1);plot(x);
title('初始语音信号时域波形');xlabel('时间轴')ylabel('幅值')
ylim([-0.2,0.2]);
subplot(2,1,2);plot(f,abs(y(1:n/2)));
title('初始语音信号频谱图');xlabel('频率Hz');ylabel('频率幅值');
ylim([0,80]);

（6）加噪语音信号的时域频域波形代码

noise=0.02*randn(1,n);
xn=x+noise';
yn=fft(xn,n);
f=fs*(0:n/2-1)/n;
figure(2)subplot(2,1,1);plot(xn);
title('语音信号加噪后的时域波形');xlabel('时间轴')ylabel('幅值')
ylim([-0.2,0.2]);
subplot(2,1,2);plot(f,abs(yn(1:n/2)));
title('语音信号加噪后的频谱图');xlabel('频率Hz');ylabel('频率幅值');
ylim([0,80]);
[x,fs]=audioread('D:\study\大三上\数字信号处理\实验\二\data1.wav');
n=length(x);
y=fft(x,n);
f=fs*(0:n/2-1)/n;

（7）设计低通滤波器并显示了滤波前后时域频域波形代码

[x,fs]=audioread( 'D:\study\大三上\数字信号处理\项目\data1.wav');
Fp=2000;Fs=1600;Ft=8000;As=100;Ap=1;
wp=2*pi*Fp/Ft;ws=2*pi*Fs/Ft;
[n,wn]=ellipord(wp,ws,Ap,As,'s');
[b,a]=ellip(n,Ap,As,wn,'s');
[B,A]=bilinear(b,a,1);
[h,w]=freqz(B,A);
% 提取音频低频分量
[x,fs]=audioread( 'D:\study\大三上\数字信号处理\项目\data1.wav');
x=x(:,1);
Y=fft(x);
y=filter(B,A,x);
Y1=fft(y);
n=0:length(x)-1;
t=(0:Fs-1)/fs;
figure(1);
subplot(2,1,1);plot(n,x);grid on;
title('低通滤波前语音信号时域波形');xlabel('频率');ylabel('幅度');
subplot(2,1,2);plot(n,y);grid on;
title('低通滤波后语音信号时域波形');xlabel('频率');ylabel('幅度');
figure(2);
subplot(2,1,1);plot(n,abs(Y));grid on;
title('低通滤波前语音信号频谱');xlabel('频率');ylabel('幅度');
axis([0 5000 0 8]);
subplot(2,1,2);plot(n,abs(Y1));grid on;
title('低通滤波后语音信号频谱');xlabel('频率');ylabel('幅度');
axis([0 5000 0 8]);

（8）设计高通滤波器并显示滤波前后时域频域波形代码

fn = 10000;%采样频率
fp = 1500;%通带截止频率
fs = 100;%阻带截止频率
Rp = 2;%通带最大衰减
Rs = 50;%阻带最小衰减
figure(1);
Wp = fp/(fn/2);
Ws = fs/(fn/2);
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[b,a] = butter(n,Wn,'high');
[H,F] = freqz(b,a,900,10000);
subplot(211);plot(F,20*log10(abs(H)))
title('高通滤波器的频率响应');xlabel('频率');ylabel('幅度');
axis([0 4000 -30 3]);
subplot(212);
pha = angle(H)*180/pi;
plot(F,pha);
title('高通滤波器的频率响应');xlabel('频率');ylabel('相位');
grid on;
[x,f]=audioread( 'D:\study\大三上\数字信号处理\项目\data1.wav');
x=x(:,1);
Y=fft(x);
y=filter(b,a,x);
Y1=fft(y);
n=0:length(x)-1;
t=(0:fs-1)/f;
figure(1);subplot(2,1,1);plot(n,x);grid on;
title('高通滤波前语音信号时域波形');xlabel('频率');ylabel('幅度');
subplot(2,1,2);plot(n,y);grid on;
title('高通滤波后语音信号时域波形');xlabel('频率');ylabel('幅度');
figure(2);
subplot(2,1,1);plot(n,abs(Y));grid on;
title('高通滤波前语音信号频谱');xlabel('频率');ylabel('幅度');
axis([0 5000 0 8]);
subplot(2,1,2);plot(n,abs(Y1));grid on;
title('高通滤波后语音信号频谱');xlabel('频率');ylabel('幅度');
axis([0 5000 0 8]);