北邮数字信号处理上机实验二答案-魔法少女雪殇

实验二：

数字信号的FFT分析数字信号处理的一个重要分支就是信号分析，而信号分析的基本工具是离散傅立叶变换。利用傅立叶变换和级数所形成的频谱分析技术作为处理连续信号的重要工具已经应用得很久了,1956年库力（Cooley）和图基（Tukey）所发展的近似频谱的快速算法为频谱分析的数字信号的谱分析铺平了道路。因此，DFT（FFT）得到广泛应用。本次实验设计了两个内容：

(1)离散信号的频谱分析

假设信号x(n)由下述信号组成：

这个信号有两根主谱线0.3pi和0.302pi靠的非常近，而另一根谱线0.45pi的幅度很小，请选择合适的长度N和窗函数，用DFT分析其频谱，得到清楚的三根谱线。要求用截图展示你的M文件，并将输入信号和分析结果信号在一张图中展示。

解题：

clc;
clear all;
N = 1000;
n = [1:1:N];
x = 0.001*cos(0.45*n*pi)+sin(0.3*n*pi)-cos(0.302*n*pi-pi/4);
y = fft(x,N);
abs = abs(y);%¶Ôfft½á¹ûÇóÄ£
w =2*pi/N*[0:1:N-1]; %Çów
subplot(2,1,1);stem(n,x);%Âö³å
subplot(2,1,2);stem(w/pi,abs);axis([0.2 0.5 0 2]);%½ÇÆµÂÊ

(2)DTMF信号频谱分析

用计算机声卡采用一段通信系统中电话双音多频（DTMF）拨号数字0～9的数据，采用快速傅立叶变换（FFT）分析这10个号码DTMF拨号时的频谱。

附录中提供了4个音频信号，请分析出这四个信号各是什么按键音，要求用截图展示你的M文件，并将四个信号在一张图中展示。

解题：

分析：

原理：滤波算法需要DTMF方案中存在频率所需的带通滤波器，可以利用Matlab滤波器设计工具fdatool创建一个具有任意截止频率的巴特沃斯带通滤波器，然后将滤波器导出至m文件

滤波器代码：

function Hd = myfilter(Fc1,Fc2)
%MYFILTER Returns a discrete-time filter object.
 
% MATLAB Code
% Generated by MATLAB(R) 9.1 and the DSP System Toolbox 9.3.
% Generated on: 09-Dec-2018 18:46:05
 
% Butterworth Bandpass filter designed using FDESIGN.BANDPASS.
 
% All frequency values are in Hz.
Fs = 8000;  % Sampling Frequency
 
N   = 10;   % Order
 
% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.bandpass('N,F3dB1,F3dB2', N, Fc1, Fc2, Fs);
Hd = design(h, 'butter');
 
% [EOF]

之后，在主文件中，为每个可能的频率创建数组。并且为每个频率创建一个滤波器，通带为频率为±25 Hz。将来自wav文件的数据作为每个滤波器的输入。每个滤波器的输出被放入一个矩阵中。为了进一步对信号进行条件化，对每个通道进行平方，然后利用均值函数对信号进行平滑。接下来，算法必须确定每个音调的起始位置与截止位置，为实现该店，做了另一个循环，他将电流值与阈值进行比对，看信号再每个点上是否工作或者停止，当他确定为一个有源脉冲时，检查滤波器组的哪两个通道由最大值，对于每一个通道有较大值的话，就进行计数，脉冲结束时，将识别计数最大的两个频率作为结果，然后用表引用结果得到字符，将该字符加入序列中，最后输出。

对于音频的噪音进行滤波，我选择的是滤波器组的方式进行滤波，它对噪声具有相对的鲁棒性，由于60Hz干扰是一个现实世界常见的干扰，所以增加了一个60hz的信号，可以有效的对白噪声进行降噪滤波。

clc;
clear all;
path = 'Z:\matlab\bin\dsp\DTMF-Decoder\';
D = dir([path,'*.wav']); %读取目录下全部wav文件
filename = cell(length(D),1);
for ii = 1:length(D)
filename = {D(ii).name};
freqs = [697,770,852,941,1209,1336,1477];
keys = ['1','2','3';'4','5','6';'7','8','9';'*','0','#'];
fs = 8000;
% 信号最大值
silence_diff_threshold = 5;
% 读取文件
[y,Fs] = audioread(char(filename));
% 减少取样
if Fs ~= fs
   disp('not match sample rate, looking for 8kHz');
   if(Fs>fs)
       disp('down')
      decimate(y,Fs/fs);
   else
       disp('quiting');
       stop;
   end
end

%figure(1);
%spectrogram(y,50,25,2048,fs,'yaxis'); %绘图
filtered_output = zeros(length(y),length(freqs));
% 循环应用于带通滤波
temp = zeros(length(y),1);
% 停止增益
Rp = 3;
Rs = 40;
for i = 1:length(freqs)
   %停止+-25Hz的距离所需频率
   Hd = myfilter(freqs(i)-25,freqs(i)+25);
   temp = filter(Hd,y);
   % remove component
   temp = temp.*temp;
   % 平均值进行平滑
   filtered_output(:,i) = movmean(temp,200);
end
subplot(4,1,ii);
plot(filtered_output,'DisplayName','filtered_output');
% 缩放不同文件的阈值
silence_threshold = max(max(filtered_output))/silence_diff_threshold;
guesses = [];%这里用来存储频率
% 七个中最大值
total1 = zeros(4,1);
total2 = zeros(3,1);
% 当其为0，用来寻找下降沿
state = 0;
%循环时间
for i = 1:length(filtered_output)
    % 根据阈值检查过滤器的总和
    if sum(filtered_output(i,:)) < silence_threshold
        % silence is found
        if state == 0
            % 这里非第一次停止
            % 重置数据，准备识别下一次音调
            total1 = zeros(4,1);
            total2 = zeros(3,1);
        else
            % 第一次停止
            % 存储上一个音调数据
            [m,f1] = max(total1);
            [m,f2] = max(total2);
            guesses = [guesses,keys(f1,f2)];
            state=0;
        end
        
    else
        % 确定过滤器哪个值最大
        % 上下波段
        [m,f1] = max(filtered_output(i,1:4));
        [m,f2] = max(filtered_output(i,5:end));
        % 增加频率总值
        total1(f1)=total1(f1)+1;
        total2(f2)=total2(f2)+1;
        % 改变状态表示有数据进行存储
        % 存储
        state = 1;
        % 检查该声音是文件最后一个
        if i==length(filtered_output)-1
            % 存储最终数据，因为此时音频可能并未结束
            % 文件结尾
            [m,f1] = max(total1);
            [m,f2] = max(total2);
            guesses = [guesses,keys(f1,f2)];
            state=0;
        end
    end
end
% 输出
sound(y);
fprintf('第%d个音频的号码识别为 .\n',ii);
disp(guesses);

end