数字信号处理实验一DFT&FFT+实验二音频处理

MATLAB 数字信号处理

一、DFT&FFT

mydft

function [Xk] = mydft(xn)
tic;
    N=length(xn);
    n = [0:1:N-1];
    k = n;
    WN= exp(-1i*2*pi/N);
    nk = n\' * k;                
    Xk = xn(1:N) * WN.^nk; 
toc;
end

myfft【这个程序有问题但是不知道哪有问题】

function [Y]=myfft(xn,N)
tic
x_odd=xn(1:2:end); %取奇数项
x_even=xn(2:2:end);%取偶数项
X1=mydft(x_even); 
X2=mydft(x_odd);
k=0;
Z=zeros(1,N);
for k=1:N/2
    Y1=X2(k)*exp(-2*pi*1i*k/N);%蝶算,奇数项提取公因子后的偶数项
    Xk1=X1(k)+Y1;
    Xk2=X1(k)-Y1;
    Z(k)=Xk1;
    Z(k+N/2)=Xk2;
end
Y=Z;
toc
end

FFT参考例程

%基2 DIT-FFT运算的MATLAB程序
function [Xk]=DITFFT(xn)
clc;close;
format compact;
%输入数据并计算常量
tic
%xn=(1:2^10);
M=nextpow2(length(xn));%用nextpow2来填充传递给fft的信号。
%这样做可以在信号长度不是2的精确幂次时加速FFT的计算,M为最靠近的幂次
N=2^M;
for m=0:N/2-1; %旋转因子指数范围
    WN(m+1)=exp(-j*2*pi/N)^m; %计算旋转因子，m为拆分后的幂次
end
A=[xn,zeros(1,N-length(xn))];%数据输入
%disp(\'输入到各存储单元的数据：\'),disp(A);
%数据倒序操作
J=0;%给倒序数赋初值
for I=0:N-1;%按序交换数据和计算倒序数
    if I<J; %条件判断及数据交换
        T=A(I+1);
        A(I+1)=A(J+1);
        A(J+1)=T;
    end
    %算下一个倒序数
    K=N/2;
    while J>=K;
        J=J-K;
        K=K/2;
    end
    J=J+K;
end
%disp(\'倒序后各存储单元的数据:\'),disp(A);
%分级按序依次进行蝶算
for L=1:M; %分级计算
   % disp(\'运算级次：\'),disp(L);
    B=2^(L-1);
    for R=0:B-1; %各级按序蝶算
        P=2^(M-L)*R;
        for K=R:2^L:N-2;
            T=A(K+1)+A(K+B+1)*WN(P+1);
            A(K+B+1)=A(K+1)-A(K+B+1)*WN(P+1);
            A(K+1)=T;
        end
    end
   % disp(\'本级运算后各存储单元的数据：\'),disp(A);
end
%disp(\'输出各存储单元的数据:\'),Xk=A;
toc
tic
%disp(\'调用FFT函数运算的结果：\')
fftxn=fft(xn,N);
toc

二、音频文件分析

function SF_plot
close;clc;clear
filepath=input(\'请输入要分析的音频文件路径:\',\'s\')
[y,fs]=audioread(filepath); 
info=audioinfo(filepath)
%sound(y,fs);
T=1/fs;
t=(0:length(y)-1)*T;
f=(0:length(y)-1)*fs/length(y);
figure(1);
yz=y(:,1);

subplot(3,1,1);
plot(t,yz);
title(\'原始信号时域\');
xlabel(\'时间\');
ylabel(\'振幅\');

subplot(3,1,2);
n=length(yz);%进行变换的点数
y1=fft(yz,n); %对n点进行傅里叶变换到频域
F=fs/length(yz); %谱分辨率，频谱间隔
plot(f,abs(y1));%左声道频谱图
set(gca,\'xticklabel\',get(gca,\'xtick\')) % 将x显示的字符替换为实际大小
title(\'原始信号频谱\');
xlabel(\'F(Hz)\');
ylabel(\'H(jw)\');
grid on

subplot(3,1,3);
n=length(yz);%进行变换的点数
y1=fft(yz,n); %对n点进行傅里叶变换到频域
F=fs/length(yz); %谱分辨率，频谱间隔

bar(f,abs(y1));%左声道频谱图
set(gca,\'xticklabel\',[0:100:1100]) % 将x显示的字符替换为实际大小
title(\'原始信号频谱\');
xlabel(\'F(Hz)\');
ylabel(\'H(jw)\');
grid on

figure(2)
subplot(211);
Pyy=y1(1:1100+1);%按照处理向量的方式得到所需声音信号的范围
f=(1:1100+1);
plot(f, 20*log10(abs(Pyy)))   %能量频谱图
xlabel(\'Frequency (Hz)\')   
ylabel(\'Power (dB)\')  
subplot(212)
bar(f,abs(Pyy)) %能量坐标图
end