Считывание исходного звукового файла

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Цель работы

- изучение программного продукта MatLab как инструмента по цифровой обработке сигналов;

- синтез цифрового фильтра с заданным АЧХ;

- ввод информации с аудио файла и его цифровая обработка (фильтрация, получение спектральной плотности с помощью ДПФ или БПФ, модуляция);

- получение спектрограмм модулированных сигналов.

Ход работы

1. Запустить программу MatLab.

Пример окна запуска версии MatLab 7.11.0.

По варианту задания произвести считывание аудио-файла в MatLab.

Для считывания, аудио-файл должен располагаться в рабочем каталоге MatLab, - например, …\ \work\test.wav, где test.wav - тестовый аудио файл (запись моно, 16 разрядов, 8 кГц, 10000 отсчетов). Установить в строке «Current folder» рабочий каталог и выполнить считывание:

>> hfile = 'test.wav';

>> y=wavread (hfile);

или одной командой

>> y=wavread ('test.wav');

возможно считывание следующими командами:

>> hfile = 'test.wav';

>> [y, Fs, bits]=wavread (hfile);

или одной командой

>> [y, Fs, bits]=wavread ('test.wav');

где Fs вернет значение частоты дискретизации, bits - количество бит для одного отсчета.

Произвести просмотр файла можно с использованием окна по команде:

>> sptool;

Для просмотра файла необходимо выполнить Файл -> Импорт… (File -> Import…)

В окне:

необходимо выбрать y; нажать на кнопку для поля «Данные» («Data») и установить 8000 для поля «Частота выборки» («Sampling Frequency»).



Затем нажать «Ок». В окне SPTool «Сигналы» («Signals») выбрать sig1[vector] и нажать «Вид» («View») под окном редактирования «Сигналы» («Signals»).

В результате на экране появится временной график сигнала, считанного из файла test.wav.

Закрыть окно графика и SPTool, не сохраняя сессию.

звуковой файл цифровой фильтр

2. По варианту задания произведем синтез полосового фильтра.

Для синтеза стандартных фильтров можно использовать:

- butter(n, w0, type) - расчет фильтров Баттерворта;

- cheby1(n, Rp, w0, type) - расчет фильтров Чебышева первого рода;

- cheby2(n, Rp, w0, type) - расчет фильтров Чебышева второго рода;

- ellip(n, Rp, Rs, w0, type) - расчет фильтров Кауэра;

где n - порядок фильтра, Rp и Rs - уровень пульсаций АЧХ в полосе пропускания и задержания (дБ), w0 и type - тип фильтра:

ФНЧ - w0 - скаляр; type - нет; (фильтр низких частот);

ФВЧ - w0 - скаляр; type - `high'; (фильтр высоких частот);

ПФ - w0 - двухэлементный вектор [w1 w2]; type - нет; (полосовой фильтр);

РФ - w0 - двухэлементный вектор [w1 w2]; type - `stop'; (режекторный фильтр).

В зависимости от количества возвращаемых параметров формируются формы цифровых фильтров:

[b, a]= - коэффициенты рекурсивной и нерекурсивной части цифрового фильтра (числителя и знаменателя функции передачи);

[z, p, k]= - нули и полюса;

[A, B, C, D]= - параметры пространства состояний.

Минимальный порядок фильтра можно вычислить, используя функции:

[n, Wn]= buttord(Wp, Ws, Rp, Rs)

[n, Wn]= cheby1ord(Wp, Ws, Rp, Rs)

[n, Wn]= cheby2ord(Wp, Ws, Rp, Rs)

[n, Wn]= ellipord(Wp, Ws, Rp, Rs)

где Wp и Ws - задают границы полос пропускания и задержания (нормированные к частоте Найквиста):

ФНЧ - Wp<Ws;

ФВЧ - Wp>Ws;

ПФ - Ws(1)< Wp(1)<Wp(2)<Ws(2);

РФ - Wp(1)< Ws(1)<Ws(2)<Wp(2).

Пусть для сигнала test.wav необходимо синтезировать ПФ с частотой среза 0,7 кГц и 2,1 кгц; с неравномерностью АЧХ 1 дБ; скоростью нарастания затухания 0,2 дБ/Гц; с затуханием в полосе непропускания 60 дБ.

Получим модель АЧХ (по затуханию), шкала частот в килогерцах и нормированная шкала:

>> Ws = [0.0625 0.2875];

>> Wp = [0.0875 0.2625];

>> Rp = 1;

>> Rs = 60; (K=0.001)

>> [n, Wn]= ellipord(Wp, Ws, Rp, Rs);



Для полученных значений параметров синтезируем фильтр командой:

>> [b, a] = ellip(n, Rp, Rs, Wn);

Для получения АЧХ и ФЧХ в логарифмической шкале можно выполнить:

>> freqz(b, a);

>> plot(w, abs(h));

Выводится график.

>> grid on;

Формируется сетка или АХЧ в линейном масштабе:

>> [h, w]=freqz(b, a);

>> plot(w, abs(h));

Выводится график.

>> grid on;

Формируется сетка или АХЧ в децибелах:

>> [h, w]=freqz(b, a);

>> plot(w, 20*log10(abs(h)));

Получим график АЧХ:

Для получения сетки на графике АЧХ, в главном окне выполним команду:

>> grid on;

Как видно с АЧХ полосовой фильтр получился с заданными параметрами.

Выполнить закрытие окна с графиком АЧХ.

Для сравнения можно получить коэффициенты рекурсивного фильтра по кусочно-линейной АЧХ строкой вида, после задания параметров для синтеза фильтра (n, f, m):

>> [b, a] = yulewalk(n, f, m);

где a и b - коэффициенты рекурсивной и нерекурсивной части цифрового фильтра; n - порядок фильтра; f - значения нормированных частот и m - значения затухания для соответствующих точек частот.

Вектора a и b будут иметь размерность n+1.

Пример синтеза фильтра по кусочно-линейной АЧХ (для получения численных значений используем соответствующие точки графика, см. выше):

>> f = [0 0.0625 0.0875 0.2625 0.2875 1];

>> m = [0.001 0.001 1 1 0.001 0.001]; (A=60 дБ - уровень затухания на графике)

>> nl = 51;

>> [b, a] = yulewalk(nl, f, m);

Для получения синтезируемого АЧХ выполним (по оси Y коэффициент передачи):

>> plot(f,m);

Для получения АХЧ в децибелах выполним:

>> [h, w]=freqz(b, a);

>> plot(w/pi, 20*log10(abs(h)));

Как видно по АЧХ синтез фильтра был выполнен успешно.

3. Выполним фильтрацию загруженного в вектор y сигнала с файла test.wav.

>> y1=filter(b, a, y);

Для просмотра результатов в визуальной форме выполним:

>> sptool;

дальнейшие действия в окне SPTool выполняем по аналогии действий пункта 1 для массива y1, вместо y. В результате получим график вида:

Как видно форма сигнала значительно изменилась.

Закрыть все активные окна, кроме основного окна программы MatLab.

4. Получим спектральную плотность сигнала для различных временных участков.

Для различных вариантов задания выбрать три различные временные точки для анализа.

Для нашего примера, выберем две начальные временные точки t1=0,12сек и t2=0,25сек, для обработки используем 512 отсчетов сигнала (кратное степени 2). Для нашего файла Fs=8000 и bits=16, тогда

для t1: N0=8000*0,12=960;

для t2: N1=8000*0,25=2000.

Выполним команды для временной точки t1:

>> nt=960;

>> Y=fft(y(nt:nt+512),512);

>> Pyy=Y.*conj(Y)/512;

>> f=8000*(1:256)/512;

>> plot(2*f, Pyy(1:256));

Выполним команды для временной точки t2:

>> nt=2000;

>> Y=fft(y(nt:nt+512),512);

>> Pyy=Y.*conj(Y)/512;

>> f=8000*(1:256)/512;

>> plot(2*f, Pyy(1:256));

Как видно из графиков спектр для точки t2 не имеет ярко выраженных спектральных «всплесков» (формант) по сравнению со спектром для точки t1.

Получим для сравнения спектры сигнала в этих точках после прохождения сигнала через полосовой фильтр (ellip), массив y1:

>> nt=960;

>> Y=fft(y1(nt:nt+512),512);

>> Pyy=Y.*conj(Y)/512;

>> f=8000*(1:256)/512;

>> plot(2*f, Pyy(1:256));

И для nt=2000:

>> nt=2000;

>> Y=fft(y1(nt:nt+512),512);

>> Pyy=Y.*conj(Y)/512;

>> f=8000*(1:256)/512;

>> plot(2*f, Pyy(1:256));



Как видно из графиков АЧХ фильтром, осуществлено выделение частот в полосе пропускания полосового фильтра, и их подавление вне пропускания полосового фильтра.

Главное окно после ранее выполненных операций имеет вид:

5. Выполним модуляцию исходного сигнала и сигнала после обработки цифровым фильтром.

При выполнении лабораторной работы модуляция сигнала производится по вариантам задания.

Для примера выберем амплитудную модуляцию (АМ), выполнив команды:

>> Fs=8000;

>> Fc=2000;

>> y2=amod(y, Fc, Fs, 'amdsb-tc');

или для версии MATLAB 7

>> y2=ammod(y, Fc, Fs);

Для просмотра результата во временной области можно аналогично пункту 1 использовать «sptool». При этом можно посмотреть и АЧХ сигнала, при условии, что частота дискретизации сигнала равна Fs/2, т.е. 4 000 Гц.

Для этого выполняем команду:

>> sptool;

В дальнейшем выполняем действия для массива y2, частота дискретизации 8000, нажимаем кнопку «Создать» («Create») под окном редактирования «Спектр» («Spectra») и кнопку «Принять» («Apply») в окне вывода спектра, что приведет к выводу окна спектра:

Шкала по оси Y на графике получена в логарифмическая (в децибелах). Для получения линейной шкалы выполним следующее:

- выберем меню «Опции» («Options»);

- подменю «Размерность Шкалы» («Magnitude Scale») и «Линейная» («Linear»), вместо «децибелы» («decbels»).

Как видно из анализа спектра после модуляции получили несущую частоту и две боковые сигнала.

Можно получить спектрограмму модулированного сигнала в виде развертки спектральной плотности по временной оси, выполнив:

>> specgram(y2, [], Fs);

Как видно из спектрограммы на частоте 2000 Гц расположена несущая, частоты выше и ниже 2000 Гц формируют боковые полосы (темно-красный цвет соответствует максимальной энергии, а темно-синий - минимальной).

Для осуществления цифровой модуляции необходимо получить многоуровневый цифровой сигнал из исходного сигнала который представляет собой вектор вещественных чисел. Количество уровней должно быть целым числом в диапазоне от 0 до М-1, где М - число уровней манипуляции (целое). Предположим, что М=8, а количество разрядов дискретизации 12. Тогда, выполнив ниже следующие команды, получим исходный массив для цифровой модуляции:

>> Bit=12;

>> Nbit = pow2(Bit)-1;

>> Amax=max(y);

>> Amin=min(y);

>> Dy=uint16(round((y-Amin)*Nbit/(Amax-Amin)));

>> N= size(Dy);

>> for k=1:1:N

yM(4*(k-1)+1)= double(bitand(Dy(k),hex2dec('0007')));

yM(4*(k-1)+2)= double(bitand(bitshift(Dy(k),-3),hex2dec('0007')));

yM(4*(k-1)+3)= double(bitand(bitshift(Dy(k),-3),hex2dec('0007')));

yM(4*(k-1)+4)= double(bitand(bitshift(Dy(k),-3),hex2dec('0007')));

end;

>>

Bit - количество разрядов дискретизации;

Nbit - два в степени 12 минус 1;

Amax и Amin - максимальные и минимальные значения исходного массива;

Dy - промежуточный без знаковый целочисленный 16-и разрядный массив (вектор);

N - количество элементов в массиве;

yM - исходный цифровой массив для модуляции, длина массива увеличивается в четыре раза по отношению к аналоговому массиву, т.к. восемь уровней (М=8) требует использования трех разрядов из 12 и тогда на один отчет аналогового сигнала приходится 12/3=4 новых отчета, что увеличивает частоту дискретизации с 8 до 32 кГц.

Для примера возьмем ЦАМ и выполним ее:

>> Fc=2000;

>> Fd=8000;

>> Fs=32000;

>> yMod=dmod(yM, Fc, Fd, Fs, 'ask', 8);

Для просмотра АЧХ и спектра сигнала можно воспользоваться «sptool» для массива yMod, задав частоту дискретизации 32000.

или получить спектрограмму модулированного сигнала, выполнив:

>> specgram(yMod, [], Fs);



Как видно из спектрограммы на частоте 2000 Гц расположена несущая, частоты выше и ниже2000 Гц формируют боковые полосы с определенной периодичностью и затуханием по частоте.

Главное окно после ранее выполненных операций имеет вид:

Выполнение работы

Для варианта задания необходимо произвести считывание исходного звукового файла и произвольно выбрать три характерные временные точки для анализа. Выбрать произвольно тип цифрового фильтра. Синтезировать цифровой фильтр и произвести фильтрацию считанного файла. По варианту выполнить модуляцию исходного и фильтрованного файла.

Содержание отчета

- цель работы;

- временной график считанного файла с указанием трех выбранных для анализа точек;

- АЧХ исходного сигнала в выбранных точках, полученных с помощью ДПФ или БПФ;

- описание выбранного цифрового фильтра и значение полученных коэффициентов функции передачи для соответствующего варианта;

- АЧХ синтезируемого фильтра;

- временной график сигнала после фильтра;

- АЧХ сигнала после фильтра в выбранных точках, полученных с помощью ДПФ или БПФ;

- спектр и спектрограмма сигнала после модуляции в одной из точек при использовании в качестве модулируемого сигнал исходный сигнал и сигнал после фильтра;

- выводы по проделанной работе.

Выводы должны содержать описание сравнения различных графиков, комментарии по спектрам и спектрограммам сигналов и т.д.

Варианты заданий

Таблица 1 - частот фильтров при скорости затухания 0,3 дБ/Гц и затухания в полосе подавления 60 дБ.

№ ПФ	1	2	3	4	5	6	7	8	9	10
кГц	0,3-1,0	0,6-1,8	0,8-2,2	1,0-2,6	1,2-2,8	1,8-3,0	2,0-3,2	2,2-3,6	2,6-3,8	3,0-4,0
Файл	U1	U2	U3	U4	U5	U1	U2	U3	U4	U5

Таблица 2 - видов модуляции.

№ мод	1	2	3	4	5	6	7	8
кГц Fc	2	2	2	2	32/8	16/8	32/8	16/8	Fs/Fd
Вид модул.	АМ DSB	АМОБ SSB	ФМ PM	ЧМ FM	ЦАМ ASK	ЦФМ PSK	ЦЧМ FSK	МЧМ MSK

Таблица вариантов.

№ вар	1	2	3	4	5	6	7	8	9	10
№ таб1	1	2	3	4	5	6	7	8	9	10
№ таб2	1	2	3	4	5	6	7	8	1	2
№ вар	11	12	13	14	15	16	17	18	19	20
№ таб1	1	2	3	4	5	6	7	8	9	10
№ таб2	3	4	5	6	7	8	1	2	3	4
№ вар	21	22	23	24	25	26	27	28	29	30
№ таб1	1	2	3	4	5	6	7	8	9	10
№ таб2	5	6	7	8	1	2	3	4	5	6

Тип фильтра выбирать произвольно.

Таблица. Функции модуляции.

Вид модуляции	Функция	Описание
АМ (DSB)	y=amod(x, Fc, Fs, `amdsb-tc', offset) или для MATLAB 7 _________________________ y=ammod(x, Fc, Fs); y = ammod(x,Fc,Fs,ini_phase) y=ammod(x,Fc,Fs,ini_phase,carramp)	amdsb-tc - тип модуляции, x - входной сигнал, Fc - несущая, Fs - частота дискретизации, offset - постоянная составляющая для моделируемого сигнала (необязательно), y - выходной сигнал. __________________________ ini_phase - начальная фаза в радианах carramp - амплитуда несущей на выходе модулятора
АМОБ (SSB)	y=amod(x, Fc, Fs, `amssb') или для MATLAB 7 _________________________ y = ssbmod(x,Fc,Fs) y = ssbmod(x,Fc,Fs,ini_phase) y=ssbmod(x,fc,fs,ini_phase,'upper')	Аналогично АМ. Используется преобразование Гильберта. _________________________ Выделяет нижнюю боковую Upper - используется для выделения верхней боковой
ФМ (PM)	y=amod(x, Fc, Fs, `pm', deviation) или для MATLAB 7 _________________________ y = pmmod(x,Fc,Fs,phasedev) y=pmmod(x,Fc,Fs,phasedev,ini_phase)	Аналогично АМ. deviation - множитель модулируемого сигнала (необязательно). _________________________ Fc - несущая (Гц), Fs - частота дискретизации (Гц) (Fs>=2Fc) phasedev -девиация фазы модулированного сигнала в радианах ini_phase - определяет начальную фазу модулированного сигнала, в радианах.
ЧМ (FM)	y=amod(x, Fc, Fs, `fm', deviation) или для MATLAB 7 _________________________ y = fmmod(x,Fc,Fs,freqdev) y= fmmod(x,Fc,Fs,freqdev,ini_phase)	Аналогично ФМ. _________________________ Fc - несущая (Гц), Fs - частота дискретизации (Гц) (Fs>=2Fc) freqdev - отклонение частоты, постоянное (Гц) модулированного сигнала. ini_phase - определяет начальную фазу модулированного сигнала, в радианах.
ЦАМ (ASK)	y=dmod(x, Fc, Fd, Fs, `ask', M) Входной сигнал х должен быть целым числом в диапазоне от 0 до М-1. или для MATLAB 7 _________________________ modem	M - число уровней манипуляции (целое), Fd - частота следования символов, Fs - частота дискретизации. Fs/Fd - положительное целое число. Fd - может быть двух элементным: первый - символьная скорость; второй - смещение (offset)/ _________________________ набор функций для реализации модуляторов и демодуляторов: DPSK, MSK, OQPSK, PSK, PAM, QAM, General QAM
ЦФМ (PSK)	y=dmod(x, Fc, Fd, Fs, `psk', M) или для MATLAB 7 _________________________ y = pskmod(x,M) y = pskmod(x,M,ini_phase) y = pskmod(x,M,ini_phase,symbol_order)	Аналогично ЦАМ.
ЦЧМ (FSK)	y=dmod(x, Fc, Fd, Fs, `fsk', M, tone) или для MATLAB 7 _________________________ y = fskmod(x,M,freq_sep,nsamp) y = fskmod(x,M,freq_sep,nsamp,Fs) y=fskmod(x,M,freq_sep,nsamp,Fs, phase_cont) y=FSKMOD(x,M,freq_sep,nsamp,Fs, phase_cont,symbol_order)	Аналогично ЦАМ. tone - шаг равномерного распределения частот вокруг Fc, так для символа 0 - f=Fc-(M-1)*tone/2, а для (M-1) - f=Fc+(M-1)*tone/2
МЧМ (MSK)	y=dmod(x, Fc, Fd, Fs, `msk')	Минимальная частотная манипуляция с разносом частот нуля и единицы на Fd/2, т.е. 0 - f=Fc-Fd/4 и 1 - f=Fc+Fd/4.

Литература

1. Цифровая обработка сигналов: учебное пособие/А.Б.Сергиенко - 3-е изд. - СПб.: БХВ-Петербург, 2011. - 768 с.: ил.

2. Основы радиоэлектроники и связи: Учебник для вузов/В.И.Нефедов - 2-е изд. Перераб. И доп. - М.: Высш. шк., 2002. - 510 с.

3. MATLAB 6/6.1/6.5 +Simulink 4/5. Основы применения. Полное руководство пользователя./Дьяконов В.П. М.: СОЛОН-Пресс. - 2002. - 768 с.

4. MATLAB 6/6.1/6.5 +Simulink 4/5 в математике и моделировании. Полное руководство пользователя./Дьяконов В.П. М.: СОЛОН-Пресс. - 2003. - 576 с.

5. MATLAB как система программирования научно-технических расчетов./В.Е.Кондрашов, С.Б.Королев. - М.: Мир, Институт стратегической стабильности Минатома РФ, 2002. - 350 с.

6. MATLAB. Анализ, идентификация и моделирование систем. Специальный справочник./Дьяконов В., Круглов В. - СПб.: Питер, 2002. - 448 с.

7. MATLAB 7.*/R2006/R2007: Самоучитель./ Дьяконов В. П. - М.: ДМК Пресс, 2008. - 768 с.: ил.

8. MATLAB 7./ Ануфриев И.Е., Смирнов А.Б., Смирнова Е.Н. - СПб.: БХВ-Петербург, 2005. - 1104 с.: ил.

9. Моделирование в системе MATLAB: Учебное пособие /Кемеровский государственный университет./ Терёхин В.В. - Новокузнецк: Куз бассвузиздат, 2004. -376с.

10. Цифровая обработка сигналов: практический подход, 2-е издание./Э.С.Айфичер, Б.У.Джервис: Пер. с англ. - М.: Издательский дом «Вильямс», 2004. - 992с.: ил.

11. Без паники! Цифровая обработка сигналов./ Юкио Сато: пер. с яп. Селиной Т.Г. - М.: Додэка-XXI, 2010. - 176 с.: ил.

Размещено на Allbest.ru

...