一、離散小波變換的音頻信号數字水印技術簡介
0 引言
近年來, 數字水印技術的作用越來越重要。數字水印技術是将一些辨別資訊直接嵌入數字載體當中, 或間接表示在信号載體中, 且不影響原載體的使用價值。通過隐藏在載體中的這些資訊, 可以判斷資訊是否被篡改, 具有防僞溯源、保護資訊安全、版權保護等作用。對于廣播轉播台站而言, 是廣播音頻的中轉站, 在廣播信号發送至千家萬戶之前務必保證信号的安全可靠, 但現在的大部分台站隻是利用人耳的判斷, 以及不同信源之間的比較, 具有較大的局限性。若利用數字水印的特性, 應用于廣播節目中可以有效地防止信号插播, 可靠地保護信号安全, 保障廣播的安全播出。
1 音頻數字水印技術分類
根據數字水印在音頻信号中的處理技術, 可将數字水印分為時域、變換域、壓縮域數字水印。
1.1 時域數字水印
在時域數字水印技術中, 直接将水印資訊嵌入至音頻信号中, 通常會選擇隐藏在信号不重要部位, 以保證其嵌入水印不影響原音頻信号的監聽效果。時域水印技術的實作較為容易且運算量小, 簡單直接, 但是魯棒性差, 容易被破解, 抵抗力較差。
1.2 變換域數字水印
在變換域數字水印中, 音頻信号需經過時域至變換域的轉換, 通常的變換域有離散餘弦變換 (DCT, Discrete Cosine Transform) 、離散傅立葉變換 (DFT, Discrete Fourier transform) 、離散小波變換 (DWT, Discrete Wavelet Transform) 等。在變換域中嵌入水印資訊, 通過反變換得到嵌入水印的音頻時域信号。變換域水印技術較時域水印技術複雜, 但變換域嵌入的水印資訊較時域而言, 不可見性更強, 隐蔽性更好, 魯棒性更好。本文的研究主要基于DWT的音頻信号的水印資訊的嵌入與提取。
1.3 壓縮域數字水印
在時域和變換域的水印技術, 都是直接将水印信号嵌入未壓縮的音頻格式中, 但是通常在音頻信号的傳輸或存儲中需要對音頻信号進行壓縮編碼 (例如WMA、MP3等) , 是以壓縮域數字水印也是水印技術也具有較大的實用價值。壓縮域數字水印技術大緻可分為三類: (1) 在非壓縮域嵌入水印, 将音頻信号與水印資訊一起壓縮; (2) 在壓縮域中, 直接将水印資訊嵌入壓縮的音頻信号中; (3) 将壓縮後的信号進行解壓縮, 然後嵌入水印資訊, 最後将水印資訊和解壓後的音頻信号一起壓縮。總的來說, 壓縮域水印技術的編解碼系統過于複雜, 受壓縮編碼格式限制大, 壓縮後的音頻信号已經去除了備援, 是以加入水印的難度大, 壓縮域水印技術有待進一步研究。
2 基于DWT的音頻水印算法
2.1 水印嵌入
本文研究的音頻水印算法是基于離散小波變換 (DWT) , 音頻信号通過DWT變換, 在變換域中嵌入水印資訊, 再經過逆變換 (IDWT) 進而得到嵌入水印的音頻信号。水印嵌入原理框圖如圖1所示。
假定水印為M1×M2的二維圖像bw, 由于音頻信号通常為一維向量, 故水印資訊在嵌入音頻信号之前需要将二維降至一維向量w, 即M=M1×M2。通常我們也可以将圖像進行打亂加密, 增強水印隐蔽性。
假定語音信号為s, 長度為N, 則s={s1, s2, s3, …, sN}由于語音信号較長在進行中一般需要進行分段, 每段長度設為N1, 故該語音信号分為K=fix (N/N1) 段進行處理, 每段語音均嵌入一個水印資訊。
小波變換是為了解決傅立葉變換的不足而提出的一種分析變換, 傅立葉變換的基函數是鋪滿整個時域的正弦信号, 對于突變信号以及變化的頻率成分資訊均不能較準确地表示。而小波變換是時間和頻率的局部變換, 更能準确地表示音頻信号的頻域特征, 常用的小波基有Haar小波、Daubechies (db N) 小波、Marr小波等。本文采用的小波基是Haar小波, 它是支撐域在t∈[0, 1]範圍内的矩形波, 定義如下:
圖1 音頻信号水印嵌入原理框圖
圖2 音頻信号水印提取原理框圖
取定Haar小波基後, 則語音信号s可以表示為:
其中Cj, k為離散小波系數, 将音頻信号分解為低頻的近似部分和高頻的細節部分, 我們在水印資訊的嵌入進行中, 主要針對代表低頻近似部分的系數向量處理, 即将水印信号放入低頻近似部分, 高頻細節部分不變, 以保證語音品質基本不變。由于嵌入的水印為二值圖像, 是以如果水印資訊的值為1, 則将對應的低頻系數增大, 相反如果值為0, 則将對應的低頻系數降低。在DWT域嵌入水印資訊後, 然後通過IDWT變換, 将語音信号變換成時域信号。
2.2 水印提取
為了保證資訊安全, 在發送端發送嵌入水印的音頻信号, 而在接收端為了确定音頻資訊的準确性, 我們通常需要提取水印以確定來源的真實性, 是以水印的提取技術也尤為重要。在水印提取過程中, 需要原始音頻信号與嵌入水印的音頻信号同時進行DWT, 再将兩者參數進行分析比較提取出水印資訊。水印提取原理框圖如圖2所示。
在前面所述的水印嵌入過程中, 将水印資訊嵌入高頻的細節部分, 是以在提取水印過程中, 我們也隻需比較原始語音信号S的低頻小波系數向量c A與嵌入水印的音頻信号s1的低頻小波系數向量c A1作比較, 若c A1>c A, 則水印資訊為1;反之則為0, 再通過向量平均, 如此得到水印資訊的一維向量, 最後通過升維得到二值圖像。
二、部分源代碼
function varargout = DSP_mark(varargin)
% DSP_MARK MATLAB code for DSP_mark.fig
% DSP_MARK, by itself, creates a new DSP_MARK or raises the existing
% singleton*.
%
% H = DSP_MARK returns the handle to a new DSP_MARK or the handle to
% the existing singleton*.
%
% DSP_MARK('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in DSP_MARK.M with the given input arguments.
%
% DSP_MARK('Property','Value',...) creates a new DSP_MARK or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before DSP_mark_OpeningFcn gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to DSP_mark_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES
% Edit the above text to modify the response to help DSP_mark
% Last Modified by GUIDE v2.5 09-Jan-2021 18:54:43
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @DSP_mark_OpeningFcn, ...
'gui_OutputFcn', @DSP_mark_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
% --- Executes just before DSP_mark is made visible.
function DSP_mark_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to DSP_mark (see VARARGIN)
% Choose default command line output for DSP_mark
handles.output = hObject;
% Update handles structure
guidata(hObject, handles);
% UIWAIT makes DSP_mark wait for user response (see UIRESUME)
% uiwait(handles.figure1);
% --- Outputs from this function are returned to the command line.
function varargout = DSP_mark_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
varargout{1} = handles.output;
% --- Executes on button press in pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDA
%--------------------------讀取原始音頻-------------------------%
cla reset;
Fs=44100;
FILE1='Rihanna - Take A Bow.wav';
[y,Fs]=audioread(FILE1);%音頻信号y，采樣率fs,采樣精度bits
y_fft=fft(y,Fs);
% y_fft_f=2*sqrt(y_fft.*conj(y_fft));
axes(handles.axes1);
plot(y);
grid on;axis tight;
title('原始音頻信号的時域波形');
xlabel('time(s)');ylabel('幅度');
axes(handles.axes2);
plot(abs(y_fft));
grid on;axis tight;
title('原始音頻信号的頻域波形');
xlabel('f(Hz)');ylabel('幅度');
sound(y,Fs);%播放原始語音
% --- Executes on button press in pushbutton2.
function pushbutton2_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
%--------------------------讀取水印音頻-------------------------%
cla reset;
Fs=44100;
FILE1='Rihanna - Take A Bow.wav';
[y,Fs]=audioread(FILE1);%音頻信号y，采樣率fs,采樣精度bits
%用db4小波對原始音頻信号進行3級小波分解
[c,l]=wavedec(y,3,'db4');%3級小波分解，低頻部分為相似，高頻部分為細節
%提取3級小波分解的低頻系數和高頻系數
ca3=appcoef(c,l,'db4',3);%提取三級小波分解的最低頻分
cd3=detcoef(c,l,3);%提取三級小波分解的次低頻分
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);
x=ca3;%提取三級小波分解的最低頻部分
%找到插入位置，檢測特征點
s=max(abs(x))*0.2;
i=find(abs(x)>s);
lx=length(x(i));
%讀取水印音頻
FILE2='test_new.wav';
[mark,Fs]=audioread(FILE2);
mark=mark(1:lx);
mark_fft=fft(mark,Fs);
% mark_fft_f=2*sqrt(mark_fft.*conj(mark_fft));
axes(handles.axes1);
plot(mark);
grid on;axis tight;
title('水印音頻信号的時域波形');
xlabel('time(s)');ylabel('幅度');
axes(handles.axes2)
plot(abs(mark_fft));
grid on;axis tight;
title('水印音頻信号的頻域波形');
xlabel('f(Hz)');ylabel('幅度');
sound(mark,Fs);
% --- Executes on button press in pushbutton3.
function pushbutton3_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton3 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
%--------------------------水印嵌入-------------------------%
cla reset;
Fs=44100;
FILE1='Rihanna - Take A Bow.wav';
[y,Fs]=audioread(FILE1);%音頻信号y，采樣率fs,采樣精度bits
%用db4小波對原始音頻信号進行3級小波分解
[c,l]=wavedec(y,3,'db4');%3級小波分解，低頻部分為相似，高頻部分為細節
%提取3級小波分解的低頻系數和高頻系數
ca3=appcoef(c,l,'db4',3);%提取三級小波分解的最低頻分
cd3=detcoef(c,l,3);%提取三級小波分解的次低頻分
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);
x=ca3;%提取三級小波分解的最低頻部分
%找到插入位置，檢測特征點
s=max(abs(x))*0.2;
i=find(abs(x)>s);
lx=length(x(i));
%讀取水印音頻
FILE2='test_new.wav';
%--------------------------讀取原始音頻-------------------------%
FILE1='Rihanna - Take A Bow.wav';
[y,Fs]=audioread(FILE1); %音頻信号y，采樣率fs,采樣精度bits
y_fft=fft(y,Fs);
% y_fft_f=2*sqrt(y_fft.*conj(y_fft));
figure(1) %繪制時域波形
plot(y);
grid on;axis tight;
title('原始音頻信号的時域波形');
xlabel('time(s)');ylabel('幅度');
figure(2) %繪制頻域波形
plot(abs(y_fft));
grid on;axis tight;
title('原始音頻信号的頻域波形');
xlabel('f(Hz)');ylabel('幅度');
sound(y,Fs); %播放原始語音
%--------------------------讀取水印音頻-------------------------%
FILE1='Rihanna - Take A Bow.wav';
[y,Fs]=audioread(FILE1); %音頻信号y，采樣率fs
%用db4小波對原始音頻信号進行3級小波分解
[c,l]=wavedec(y,3,'db4'); %3級小波分解，低頻部分為相似，高頻部分為細節
%提取3級小波分解的低頻系數和高頻系數
%這個會傳回信号y在level 3（小波分解級數）的小波分解系數
ca3=appcoef(c,l,'db4',3); %提取三級小波分解的最低頻分
cd3=detcoef(c,l,3); %提取三級小波分解的次低頻分量
cd2=detcoef(c,l,2);
cd1=detcoef(c,l,1);
x=ca3; %提取三級小波分解的最低頻部分
%找到插入位置，檢測特征點
s=max(abs(x))*0.2;
i=find(abs(x)>s);
lx=length(x(i));
%讀取水印音頻
FILE2='test_new.wav';
[mark,Fs]=audioread(FILE2);
mark=mark(1:lx);
mark_fft=fft(mark,Fs);
% mark_fft_f=2*sqrt(mark_fft.*conj(mark_fft));
figure(3)
plot(mark);
grid on;axis tight;
title('水印音頻信号的時域波形');
xlabel('time(s)');ylabel('幅度');
figure(4)
plot(abs(mark_fft));
grid on;axis tight;
title('水印音頻信号的頻域波形');
xlabel('f(Hz)');ylabel('幅度');
sound(mark,Fs); 三、運作結果
四、matlab版本及參考文獻
1 matlab版本
2014a
2 參考文獻
[1]韓紀慶,張磊,鄭鐵然.語音信号處理（第3版）[M].清華大學出版社，2019.