clc
clear all
close all
dbstop if error
%% 频率值与音阶音调标注对应，12个音阶，每个音阶有9个音调
[dataStruct, textdata, colheaders] = importdata('音调与频率.txt');
musicFreqMatrix = dataStruct.data(:, 2:end);
musicFreqNameMatrix = dataStruct.textdata(2:end);
[row, col] = size(musicFreqMatrix); % 第一列是9个级别
n = 1;
musicFreq = zeros(row*col, 1);
musicFreqName = cell(row*col, 1);
for i = 1 : row
for j = 1 : col
musicFreq(n) = musicFreqMatrix(i, j);
musicFreqName{n} = [musicFreqNameMatrix{j}, num2str(i-1)];
n = n + 1;
end
end
%% 音乐处理
[musicData, Fs]=audioread('123我爱你.mp3');
% sound(musicData, Fs);
musicData(musicData(:,1) == 0, :) = []; % 删除音乐前后空白
musicLength = length(musicData(:, 1)); % 音乐总时长
eachTime = 1 * Fs; % 数据分析使用的数据长度
% eachTime = musicLength
Count = floor(musicLength / eachTime);
%% FIR滤波，低通，滤波抽头系数计算
% N = 380; % Order
% Fpass = 800; % Passband Frequency
% Fstop = 1000; % Stopband Frequency
% Wpass = 1; % Passband Weight
% Wstop = 80; % Stopband Weight
% dens = 381; % Density Factor
%
% % Calculate the coefficients using the FIRPM function.
% b = firpm(N, [0 Fpass Fstop Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop], ...
% {dens});
% Hd = dfilt.dffir(b);
Fpass = 8000; % Passband Frequency
Fstop = 12000; % Stopband Frequency
Dpass = 0.057501127785; % Passband Ripple
Dstop = 0.0001; % Stopband Attenuation
dens = 20; % Density Factor
% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fpass, Fstop]/(Fs/2), [1 0], [Dpass, Dstop]);
% Calculate the coefficients using the FIRPM function.
b = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);
%% FIR滤波，带通（500~2000Hz）,过渡带（200Hz），阻带抑制50dB
% N = 250; % Order
% Fstop1 = 300; % First Stopband Frequency
% Fpass1 = 500; % First Passband Frequency
% Fpass2 = 2000; % Second Passband Frequency
% Fstop2 = 2200; % Second Stopband Frequency
% Wstop1 = 50; % First Stopband Weight
% Wpass = 1; % Passband Weight
% Wstop2 = 50; % Second Stopband Weight
% dens = 251; % Density Factor
% Hd = filter_FIR_bandpass(Fs, N, Fstop1, Fpass1, Fpass2, Fstop2, Wstop1, Wpass, Wstop2, dens);
figure(1)
for i = 0 : Count-1
musicDataTmp = musicData(1 + eachTime*i : eachTime + eachTime*i); % 获取录音数据
%% FIR滤波
mySpeech_FIR_filter = fftfilt(Hd.Numerator, musicDataTmp);
%% 时域转频域处理
fftN = 2^(ceil(log2(length(mySpeech_FIR_filter)))); % fft处理点数调整为2的整数次幂
sigL = length(mySpeech_FIR_filter);
% fft后直流分量部分增大了N倍、非直流分量部分增大了N/2倍，非直流分量数据点归一化需要除以N/2。
musicData_fft = fft(musicDataTmp, fftN) * 2 / sigL;
musicData_fft = musicData_fft(1 : fftN/2);
musicData_fft(1) = musicData_fft(1)/2; % 直流分量
musicData_FIR_fft = fft(mySpeech_FIR_filter, fftN) * 2 / sigL;
musicData_FIR_fft = musicData_FIR_fft(1 : fftN/2);
musicData_FIR_fft(1) = musicData_FIR_fft(1)/2; % 直流分量
%% 绘制时域波形图
subplot(2,1,1)
time = (1 : length(musicDataTmp)) ./ Fs * 1000; % ms
plot(time, musicDataTmp, 'b', time, mySpeech_FIR_filter, 'r')
legend('原始信号', 'FIR滤波')
xlabel('时间/ms')
ylabel('amp')
grid on
drawnow
%% 绘制频域图
h = subplot(2,1,2);
set(h, 'XTickLabelMode', 'manual')
set(h, 'XTickMode', 'manual')
freq = (0 : Fs/fftN : Fs/2 - Fs/fftN)'; % 生成真实频率横坐标
index = freq < 700;
plot(freq(index), db(abs(musicData_fft(index))), 'b', freq(index), db(abs(musicData_FIR_fft(index))), 'r')
ylim([-120, 0])
set(h, 'XTick', musicFreq)
set(h, 'XTickLabel', musicFreqName)
legend('原始信号', 'FIR滤波')
xlabel('音阶')
ylabel('dB')
grid on
drawnow
pause(1)
end
