audio_wm/embed.cpp

#include <iostream>
#include <vector>
#include <opencv2/opencv.hpp>
#include <filesystem>
#include <sndfile.h>
#include <fmt/core.h>
#include <windows.h>
#include <fstream>
#include <cmath>
#include <numeric>
#include <sstream>
#include <cstring>
#include <Eigen/Dense>
#include "global.h"
#include "dr_wav.h"
#include "embed.h"
#include "util.h"
#include "wavelet.h"
//#include "miniaudio.h"
#include "dr_mp3.h"

using namespace std;
using namespace cv;
namespace fs = std::filesystem;
using namespace Eigen;

const int						k1 = 5;
const double					PI = 3.1416;
int                             G_L = 0;
const int                       G_NMAX_32 = 40;
const int                       G_SIZE_32 = 2 * G_NMAX_32 + 1;
const int                       G_TOTAL_ELEMENTS_32 = G_SIZE_32 * G_SIZE_32;
vector<vector<int>>             G_ZMLIST_32(G_TOTAL_ELEMENTS_32, vector<int>(2, 0));
vector<complex<double>>         G_A_NM_32(G_TOTAL_ELEMENTS_32, (0, 0));
vector<vector<int>>             G_ZMLIST_SELECTED_32(G_TOTAL_ELEMENTS_32, vector<int>(2, 0));
const int                       G_NMAX_10 = 20;
const int                       G_SIZE_10 = 2 * G_NMAX_10 + 1;
const int                       G_TOTAL_ELEMENTS_10 = G_SIZE_10 * G_SIZE_10;
vector<vector<int>>             G_ZMLIST_10(G_TOTAL_ELEMENTS_10, vector<int>(2, 0));
vector<complex<double>>         G_A_NM_10(G_TOTAL_ELEMENTS_10, (0, 0));
vector<vector<int>>             G_ZMLIST_SELECTED_10(G_TOTAL_ELEMENTS_10, vector<int>(2, 0));
vector<complex<double>>         G_A_NM_SELECTED;
vector<complex<double>>         G_A_NM_MODIFIED;
vector<vector<complex<double>>> G_GP_REC_BEFOREMODIFY;
vector<vector<complex<double>>> G_GP_REC;

int wm_embed(const wchar_t* path, const wchar_t* save_path, const int* wm, const int wm_size, const double start_time, const double end_time) {
	if (!canFormSquareMatrix(wm_size)) return 2;
	int is_mp3 = 0;
	int Fs = 0;
	int wavelength = 0;
	//double D = 0.8; // init 0.2
	double DD = 0.0012; // init DD 0.0012
	double EE = 0.00002; // init EE 0.00002
	int begin = 4095;
	int block_length = 0;
	vector<int> syn = {1, 1, 1, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 1, 0};
	const int Lsyn = syn.size();
	fs::path fname = path;
	fs::path s_path = save_path;
	int channels;
	//int bitrate;
	vector<double> yo;
	vector<double> yo_origin;
	int start = 0;
	string audio_name = ConvertUTF16ToMBCS(path);
	string save_name = ConvertUTF16ToMBCS(save_path);

	cout << "\n read file:" << audio_name << endl;

	if (isMp3File(fname, audio_name))
	{
		// is mp3 
		is_mp3 = 1;
		long Fs_long = 0;
		yo_origin = mp3ToDoubleVector(audio_name.c_str(), channels, Fs_long, wavelength);
		Fs = (int)Fs_long;
	}
	else {
		drwav wav;
		if (!drwav_init_file(&wav, audio_name.c_str(), nullptr)) {
			std::cerr << "Error opening WAV file." << std::endl;
			return 1;
		}

		yo_origin = readWav(wav, wavelength, channels, Fs);
	}

	// stage 1 end
	if ((end_time - start_time) < 25.0)
		// The difference between start_time and end_time is not less than 25 seconds
		return 4;

	//int r = static_cast<double>(wavelength) / Fs;
	//double embed_range = ceil(r * 0.58);

	//if ((end_time - start_time) < embed_range) {
	//	cout << "[embed error] Please increase the embedding range to " << embed_range << " seconds or more." << endl;
	//	return 6;
	//}

	if (Fs < 44100) {
		wstring resample_path = AddSuffixToPath(path, L"_resample");
		vector<double> out;
		convert_wav(yo_origin, Fs, channels, resample_path.c_str());
	}

	int max_wavelength = Fs * static_cast<int>(end_time - start_time) * channels;
	// stage 1 start
	//const int frame_size = max_wavelength; // 每帧大小
	//const int hop_size = Fs; // 帧移
	//size_t num_frames = (yo_origin.size() - frame_size) / hop_size + 1;
	//std::vector<double> energy(num_frames);

	//// Calculate short-term energy
	//for (size_t i = 0; i < num_frames; ++i) {
	//	double frame_energy = 0.0;
	//	for (int j = 0; j < frame_size; ++j) {
	//		frame_energy += yo_origin[i * hop_size + j] * yo_origin[i * hop_size + j];
	//	}
	//	energy[i] = frame_energy;
	//}

	//// Calculate the mean and standard deviation of energy
	//double mean_energy = std::accumulate(energy.begin(), energy.end(), 0.0) / num_frames;
	//double sq_sum = std::inner_product(energy.begin(), energy.end(), energy.begin(), 0.0);
	//double stddev_energy = std::sqrt(sq_sum / num_frames - mean_energy * mean_energy);

	//// Set a threshold, such as adding twice the standard deviation to the mean
	//double threshold = mean_energy + 2 * stddev_energy;

	//// Detect noisy areas
	//for (size_t i = 0; i < num_frames; ++i) {
	//	if (energy[i] > threshold) {
	//		start = i * hop_size;
	//		std::cout << "Noisy frame detected at frame index: " << i << ", yo_origin[" << i * hop_size << "], sec: " << static_cast<double>(i * hop_size) / (Fs * 2) << "s,Energy: " << energy[i] << std::endl;

	//		// 将噪杂的帧添加到yo中
	//		for (int j = 0; j < max_wavelength; ++j) {
	//			yo.push_back(yo_origin[start + j]);
	//		}

	//		wavelength = max_wavelength / channels;
	//		break;
	//	}
	//}

	bool split = yo_origin.size() > max_wavelength;

	if (yo.empty()) {
		// 如果音频长度大于25s 则根据start_time, end_time 来确定要嵌入水印的范围
		// getAudio20sFront(split, max_wavelength, Fs, channels, wavelength, yo_origin, yo);

		if (split) {
			start = static_cast<int>(start_time) * Fs * channels;
			cout << "start: " << start << endl;
			splitAudioByTime(max_wavelength, Fs, channels, wavelength, yo_origin, yo, start_time, end_time, start);
		}
		else {
			yo = yo_origin;
		}
	}

	vector<vector<double>> audio_data(yo.size() / channels, vector<double>(channels, 0.0));
	cout << "yo.size():" << yo.size() << endl;
	cout << "audio_data.size():" << audio_data.size() << endl;
	if (channels == 2)
	{
		for (int i = 0; i < yo.size(); i++)
		{
			if (i % 2 == 0) {
				audio_data[i / 2][0] = yo[i];
			}
			if (i % 2 == 1) {
				audio_data[i / 2][1] = yo[i];
			}
		}
	}
	else
	{
		for (int i = 0; i < yo.size(); i++)
		{
			audio_data[i][0] = yo[i];
		}
	}
	std::vector<int> wm_vector(wm, wm + wm_size);
	//vector<vector<int>> ww_matrix = matToVector(ww, rows, cols);
	int row = static_cast<int>(std::sqrt(wm_size));
	vector<vector<int>> ww_matrix = int_ott(wm_vector, row);
	int rows = ww_matrix.size();
	int cols = ww_matrix[0].size();
	cout << "ww_rows:" << ww_matrix.size() << "\t ww_cols:" << ww_matrix[0].size() << endl;
	// 图像置乱
	vector<vector<int>> w = general(ww_matrix, rows, cols, 5, 6);
	// 降维
	vector<int> w1 = int_tto(w, rows);
	//block_length = rows * cols * 8 * 8;
	block_length = 65536;
	int i = 1;
	//int left = i * (Lsyn * k1 + block_length);
	int right = wavelength - (Lsyn * k1 + block_length) * 2;
	vector<int> t0;

	while (i * (Lsyn * k1 + block_length) < wavelength - ((Lsyn * k1 + block_length)))
	{
		int bb = begin + (i - 1) * (Lsyn * k1 + block_length);
		cout << "\n embed barker code form " << bb << " to " << bb + Lsyn * k1 << endl;

		for (int mm = 0; mm < Lsyn; mm++)
		{
			//vector<double> temp_vec_subsyn;
			int front = bb + (mm * k1);
			int back = bb + ((mm + 1) * k1);
			//计算指定范围均值
			double tempmean = calculateMean(audio_data, front + 1, back, 0); // front = 4096 back = 4100
			int temp = static_cast<int>(floor(tempmean / D));
			double tempmeanl;

			if ((temp % 2 + 2) % 2 == syn[mm])
			{
				tempmeanl = temp * D + D / 2.0;
			}
			else {
				tempmeanl = temp * D - D / 2.0;
			}
			// 更新指定范围的yo
			for (int j = front; j <= back; j++)
			{
				audio_data[j][0] += tempmeanl - tempmean;
			}
		}

		t0.push_back(begin + (i * Lsyn * k1) + ((i - 1) * block_length) + 1);
		vector<double> BLOCK1;

		for (int j = t0[i - 1]; j < t0[i - 1] + block_length; j++)
		{
			BLOCK1.push_back(audio_data[j][0]);
		}

		vector<double> result = wavelet(BLOCK1, 1);
		result.resize(65536);
		BLOCK1.clear();
		vector<vector<double>> I = ott(result, 256);
		int n = I.size();  // 行数
		int m = I[0].size();  // 列数
		// 极坐标转换相关参数
		int M = 4 * n;

		//vector<vector<double>> fp(M, vector<double>(M, 0.0));
		//vector<vector<double>> Gp(M, vector<double>(M, 0.0));
		vector<vector<double>> rr(M, vector<double>(M, 0.0));
		vector<vector<double>> oo(M, vector<double>(M, 0.0));
		// 声明和初始化 vv 和 uu 矩阵
		vector<vector<int>> vv(M, vector<int>(M, 0));
		vector<vector<int>> uu(M, vector<int>(M, 0));

		MatrixXd fp(M, M);
		MatrixXd Gp(M, M);

		// 填充 fp 和 Gp 数组
		for (int u = 0; u < M; u++) {
			for (int v = 0; v < M; v++) {
				double rrr = static_cast<double> (u) / M;
				double ooo = (2 * PI * v) / M;
				int kk = ceil(rrr * (n / 2) * sin(ooo));
				int ll = ceil(rrr * (n / 2) * cos(ooo));
				fp(u, v) = I[(-1) * kk + (n / 2)][ll + (n / 2) - 1];
				Gp(u, v) = fp(u, v) * std::sqrt(static_cast<double> (u) / (2 * M));
			}
		}
		// 填充 vv 和 uu 矩阵
		for (int k = 0; k < M; k++) {
			for (int j = 0; j < M; j++) {
				vv[k][j] = j + 1;
				uu[k][j] = k + 1;
			}
		}
		// 计算 rr 和 oo 矩阵
		for (int k = 0; k < M; k++) {
			for (int j = 0; j < M; j++) {
				rr[k][j] = static_cast<double>(uu[k][j] - 1) / M;
				oo[k][j] = static_cast<double>(2 * PI * (vv[k][j] - 1)) / M;
			}
		}

		vector<vector<double>> Gp_watermarked(M, vector<double>(M, 0.0));

		PQIMfft_embed(Gp, fp, w1, M, DD, EE, oo, rr, Gp_watermarked);

		cout << "初始化 I_incircle 为全零矩阵" << endl;

		vector<vector<double>> I_incircle(n, vector<double>(m, 0.0));
		// 嵌入水印
		cout << "\n 嵌入水印" << endl;

		for (int u = 0; u < M; u++) {
			for (int v = 0; v < M; v++) {
				double r = static_cast<double>(u) / M;
				double o = (2 * PI * v) / M;
				int kkk = static_cast<int>(std::ceil(r * (n / 2) * std::sin(o)));
				int lll = static_cast<int>(std::ceil(r * (n / 2) * std::cos(o)));
				I_incircle[-kkk + (n / 2)][lll + (n / 2) - 1] = Gp_watermarked[u][v];
			}
		}
		// 重构水印
		cout << "重构水印" << endl;
		// vector<vector<double>> I_watermarked = I_incircle;
		cout << "I_incircle's rows:" << I_incircle.size() << "I_incircle's cols: " << I_incircle[0].size() << endl;

		for (int ii = 0; ii < n; ii++) {
			for (int jj = 0; jj < m; jj++) {
				if (I_incircle[ii][jj] == 0) {
					I_incircle[ii][jj] = I[ii][jj];
				}
			}
		}

		cout << "\n go to tto;" << endl;

		vector<double> BLOCK2 = tto(I_incircle, 256);

		vector<double> new_BLOCK2(BLOCK2.size(), 0.0);

		BLOCK2.insert(BLOCK2.end(), new_BLOCK2.begin(), new_BLOCK2.end());

		cout << "\n go to inverseWavelet;" << endl;

		vector<double> signal = inverseWavelet(BLOCK2, 1);

		//cout << "\n NO." << i << " , form " << t0[i - 1] << " to " << t0[i - 1] + block_length - 1 << endl;

		for (int j = t0[i - 1]; j < t0[i - 1] + block_length; j++) {
			audio_data[j][0] = signal[j - t0[i - 1]];
		}
		
		double sub_progress = double(i * (Lsyn * k1 + block_length)) / (wavelength - (Lsyn * k1 + block_length) * 2) * 100;

		g_audio_wm_progress = sub_progress;

		cout << "\n 总完成进度:" << sub_progress << endl;
		cout << "==============================================================" << endl;

		i++;
		//left = i * (Lsyn * k1 + block_length);
	}
	cout << "t0:" << endl;
	printIntSignal(t0, 0);
	//trans audio_data to yo
	yo.clear();
	if (channels == 2) {
		for (int j = 0; j < audio_data.size(); j++)
		{
			yo.push_back(audio_data[j][0]);
			yo.push_back(audio_data[j][1]);
		}
	}
	else {
		for (int j = 0; j < audio_data.size(); j++)
		{
			yo.push_back(audio_data[j][0]);
		}
	}
	cout << "yo.size: " << yo.size() << endl;
	printSignal(yo, 10);
	cout << "start: " << start << endl;
	if (split)
	{
		for (size_t i = 0; i < yo.size(); ++i) {
			yo_origin[start + i] = yo[i];
		}
	}
	else 
	{
		yo_origin = yo;
	}
	// save audio
	//if (s_path.extension() != ".wav")
	//	s_path.replace_extension(".wav");
	const wchar_t* wchar_save_path = s_path.c_str(); // 使用完整的保存路径
	if (is_mp3) {
		save_audio_drmp3(wchar_save_path, yo_origin, Fs, channels);
	}
	else {
		save_audio_drwav(wchar_save_path, yo_origin, Fs, channels);
	}
	//fs::remove(fname);
	return 0;
}

void saveAudio(const char* outputFilename, vector<double> &yo, int sampleRate, int channels) {
	SF_INFO sfinfo{};
	sfinfo.samplerate = sampleRate;   // 设置采样率
	sfinfo.channels = channels;         // 设置声道数
	sfinfo.format = SF_FORMAT_WAV | SF_FORMAT_PCM_16;  // 设置文件格式
	// todo save audio
	SNDFILE* outfile = sf_open(outputFilename, SFM_WRITE, &sfinfo);

	sf_count_t count = sf_write_double(outfile, yo.data(), yo.size());
	cout << "\n保存地址:" << outputFilename << endl;
	if (count != yo.size()) {
		std::cerr << "Error writing frames to file" << std::endl;
		sf_close(outfile);
		return;
	}
	// 关闭文件
	sf_close(outfile);
	std::cout << "\nFile saved successfully." << std::endl;
}

vector<double> inverseWavelet(vector<double> swt_output, int level) {
	vector<double> iswt_output;

	iswt(swt_output, level, "haar", iswt_output);

	return iswt_output;
}

void PQIMfft_embed(MatrixXd Gp, MatrixXd fp,
					vector<int> &w1, int M, double DD, double EE, vector<vector<double>> &oo,
					vector<vector<double>> &rr, vector<vector<double>> &Gp_watermarked)
{
	G_L = w1.size();
	int row = static_cast<int>(std::sqrt(G_L));
	vector<vector<int>> zmlist;
	vector<vector<int>> zmlist_selected;
	vector<complex<double>> a_nm;
	int nmax = 0;

	if (row == 10) {
		nmax = G_NMAX_10;
		a_nm = G_A_NM_10;
		zmlist = G_ZMLIST_10;
	}
	else {
		nmax = G_NMAX_32;
		a_nm = G_A_NM_32;
		zmlist = G_ZMLIST_32;
	}
	zhishujufenjie62(Gp, M, nmax, a_nm, zmlist);
	// 对应的n和m
	vector<int> index_n(zmlist.size(), 0);
	vector<int> index_m(zmlist.size(), 0);

	for (int i = 0; i < zmlist.size(); i++) {
		index_n[i] = zmlist[i][0];
		index_m[i] = zmlist[i][1];
	}
	// 找到合适的嵌入位置
	vector<int> index_suitable;

	for (int i = 0; i < zmlist.size(); i++) {
		if (row == 10) {
			// 10*10 embed range
			if (((index_m[i] == 0) && (index_n[i] < 0)) ||
				((index_m[i] > 0) && (index_m[i] <= 10) && (-4 <= index_n[i]) && (index_n[i] <= 4))) {
				index_suitable.push_back(i);
			}
		}
		else {
			// 32*32 embed range
			if (((index_m[i] == 0) && (index_n[i] < 0)) ||
				((index_m[i] > 0) && (index_m[i] <= 30) && (-16 <= index_n[i]) && (index_n[i] <= 16))) {
				index_suitable.push_back(i);
			}
		}
	}

	index_n.clear();
	index_m.clear();
	// 随机选出嵌入水印的位置
	vector<int> index_selected_embed(index_suitable.begin(), index_suitable.begin() + G_L);

	index_suitable.clear();
	// 得到在zmlist中的对应的数据，为下一步重构准备的
	//G_ZMLIST_SELECTED_10.resize(G_L, vector<int>(2, 0));
	zmlist_selected.resize(G_L, vector<int>(2, 0));
	G_A_NM_SELECTED.resize(G_L, (0,0));

	for (int i = 0; i < G_L; i++) {
		//G_ZMLIST_SELECTED_10[i] = G_ZMLIST_10[index_selected_embed[i]];
		zmlist_selected[i] = zmlist[index_selected_embed[i]];
		G_A_NM_SELECTED[i] = a_nm[index_selected_embed[i]];
	}

	index_selected_embed.clear();

	G_A_NM_MODIFIED = G_A_NM_SELECTED;

	vector<double> III(G_L, 0.0);

	vector<double> IIII(G_L, 0.0);

	for (int i = 0; i < G_L; i++) {
		complex<double> A_selected = G_A_NM_MODIFIED[i];

		III[i] = ((-DD) + sqrt(pow(DD, 2) + 2 * EE * abs(A_selected))) / EE;

		double flag = 2 * floor((III[i] + 1) / 2);

		if (flag >= III[i]) {
			IIII[i] = 2 * floor((III[i] + 1) / 2) - w1[i];
		}
		else {
			IIII[i] = 2 * floor((III[i] + 1) / 2) + w1[i];
		}

		double aa = DD * IIII[i];

		double bb = EE * ((pow(IIII[i], 2)) * 0.5);

		G_A_NM_MODIFIED[i] = complex_sign(A_selected) * (aa + bb);
	}

	III.clear();
	IIII.clear();

	G_GP_REC_BEFOREMODIFY.resize(M, vector<complex<double>>(M, complex<double>(0, 0)));

	//cout << "\n go to zhishujurec622" << endl;

	//zhishujurec622(G_A_NM_SELECTED, G_ZMLIST_SELECTED_10, M, oo, rr, G_GP_REC_BEFOREMODIFY);
	zhishujurec622(G_A_NM_SELECTED, zmlist_selected, M, oo, rr, G_GP_REC_BEFOREMODIFY);

	G_A_NM_SELECTED.clear();

	//cout << "--------------------------" << endl;

	G_GP_REC.resize(M, vector<complex<double>>(M, complex<double>(0, 0)));

	//zhishujurec622(G_A_NM_MODIFIED, G_ZMLIST_SELECTED_10, M, oo, rr, G_GP_REC);
	zhishujurec622(G_A_NM_MODIFIED, zmlist_selected, M, oo, rr, G_GP_REC);

	G_A_NM_MODIFIED.clear();
	G_ZMLIST_SELECTED_10.clear();
	//水印图像
	for (int i = 0; i < M; i++) {
		for (int j = 0; j < M; j++) {
			Gp_watermarked[i][j] = fp(i, j) - std::real(G_GP_REC_BEFOREMODIFY[i][j]) + std::real(G_GP_REC[i][j]);
		}
	}

	G_GP_REC_BEFOREMODIFY.clear();
	G_GP_REC_BEFOREMODIFY.clear();
	G_GP_REC.clear();

	cout << "rows: " << Gp_watermarked.size() << "; cols: " << Gp_watermarked[0].size() << endl;
}

void zhishujurec622(vector<complex<double>> A_nm, vector<vector<int>>& zmlist_selected, int M,
	vector<vector<double>> oo, vector<vector<double>> rr, vector<vector<complex<double>>> &result) {
	// Initialize matrices
	//const int SIZE = static_cast<int>(std::sqrt(G_L));
	//int offset = 20;
	//if (SIZE == 32)
	//	offset = 40;
	//const int COUNT_SIZE = 57;
	//vector<int> count(COUNT_SIZE, 0);
	//vector<double> temp_oo_row = oo[0];
	//vector<double> temp_rr_col;
	//vector<complex<double>> RBF1_col(M, complex<double>(0.0, 0.0));
	//vector<complex<double>> RBF2_col(M, complex<double>(0.0, 0.0));
	//// Iterate through the 32x32 grid
	//for (int x = 0; x < SIZE; x++) {
	//	//start = clock();
	//	for (int y = 0; y < SIZE; y++) {
	//		int row = x * SIZE + y;
	//		int nj = zmlist_selected[row][0];
	//		int mc = zmlist_selected[row][1];
	//		complex<double> A_nm_value = A_nm[row];
	//		complex<double> conj_A_nm_value = conj(A_nm_value);
	//		int cnt_nj = count[nj + offset]; // 32*32 +40; 10*10 +20; 16*16 +20;
	//		double coefficient = 2 * nj * PI;
	//		for (int i = 0; i < M; i++)
	//		{
	//			double rr_val = rr[i][0]; // 行相同
	//			double image_rr_val = coefficient * rr_val;
	//			double cos_image_rr_val = cos(image_rr_val);
	//			double sin_image_rr_val = sin(image_rr_val);
	//			if (cnt_nj == 0)
	//			{
	//				double sqrt_part = sqrt(2.0 / rr_val);
	//				RBF1_col[i] = sqrt_part * complex<double>(cos_image_rr_val, sin_image_rr_val);
	//				RBF2_col[i] = sqrt_part * complex<double>(cos_image_rr_val, -sin_image_rr_val);
	//			}
	//			complex<double> coefficient1 = A_nm_value * RBF1_col[i];
	//			complex<double> coefficient2 = conj_A_nm_value * RBF2_col[i];
	//			for (int j = 0; j < M; j++)
	//			{
	//				double oo_val = temp_oo_row[j];// 列相同
	//				double image_rr_val = mc * oo_val;
	//				double cos_oo_val = cos(image_rr_val);
	//				double sin_oo_val = sin(image_rr_val);
	//				complex<double> front = coefficient1 * complex<double>(cos_oo_val, sin_oo_val);
	//				complex<double> behind = coefficient2 * complex<double>(cos_oo_val, -sin_oo_val);
	//				result[i][j] += front + behind;
	//			}
	//		}
	//		if (cnt_nj == 0)
	//		{
	//			count[nj + offset] = 1; // 32*32 +40; 10*10 +20; 16*16 +20;
	//		}
	//	}
	//}
	
	// ========================================================================================================================
	// Initialize matrices
	const int SIZE = static_cast<int>(std::sqrt(G_L));
	int offset = 20;
	if (SIZE == 32)
		offset = 40;
	const int COUNT_SIZE = 57;
	MatrixXd eigen_oo(oo.size(), oo[0].size());
	for (size_t i = 0; i < oo.size(); ++i) {
		for (size_t j = 0; j < oo[i].size(); ++j) {
			eigen_oo(i, j) = oo[i][j];
		}
	}

	MatrixXcd eigen_A_nm(A_nm.size(), 1);
	for (size_t i = 0; i < A_nm.size(); ++i) {
		eigen_A_nm(i, 0) = A_nm[i];
	}

	MatrixXd eigen_rr(rr.size(), rr[0].size());
	for (size_t i = 0; i < rr.size(); ++i) {
		for (size_t j = 0; j < rr[i].size(); ++j) {
			eigen_rr(i, j) = rr[i][j];
		}
	}
	VectorXi count = VectorXi::Zero(COUNT_SIZE);
	VectorXd temp_oo_row = eigen_oo.row(0);
	VectorXcd RBF1_col(M);
	VectorXcd RBF2_col(M);
	// Iterate through the 32x32 grid
	for (int x = 0; x < SIZE; x++) {
		//start = clock();
		for (int y = 0; y < SIZE; y++) {
			int row = x * SIZE + y;
			int nj = zmlist_selected[row][0];
			int mc = zmlist_selected[row][1];
			// complex<double> A_nm_value = A_nm[row];
			complex<double> A_nm_value = eigen_A_nm(row);
			complex<double> conj_A_nm_value = conj(A_nm_value);
			int cnt_nj = count(nj + offset);
			double coefficient = 2 * nj * PI;
			for (int i = 0; i < M; i++)
			{
				double rr_val = eigen_rr(i, 0);
				double image_rr_val = coefficient * rr_val;
				double cos_image_rr_val = cos(image_rr_val);
				double sin_image_rr_val = sin(image_rr_val);
				if (cnt_nj == 0)
				{
					double sqrt_part = sqrt(2.0 / rr_val);
					RBF1_col(i) = sqrt_part * complex<double>(cos_image_rr_val, sin_image_rr_val);
					RBF2_col(i) = sqrt_part * complex<double>(cos_image_rr_val, -sin_image_rr_val);
				}
				complex<double> coefficient1 = A_nm_value * RBF1_col[i];
				complex<double> coefficient2 = conj_A_nm_value * RBF2_col[i];
				for (int j = 0; j < M; j++)
				{
					double image_rr_val = mc * temp_oo_row(j);
					double cos_oo_val = cos(image_rr_val);
					double sin_oo_val = sin(image_rr_val);
					complex<double> front = coefficient1 * complex<double>(cos_oo_val, sin_oo_val);
					complex<double> behind = coefficient2 * complex<double>(cos_oo_val, -sin_oo_val);
					result[i][j] += front + behind;
					//result[i][j] += coefficient1 * complex<double>(cos(mc * temp_oo_row(j)), sin(mc * temp_oo_row(j))) +
					//	coefficient2 * complex<double>(cos(mc * temp_oo_row(j)), -sin(mc * temp_oo_row(j)));
				}
			}
			if (cnt_nj == 0)
			{
				count(nj + offset) = 1;
			}
		}
	}
}

// 对于复数的sign函数，返回其实部和虚部的符号组成的复数  
template<typename T>

std::complex<T> complex_sign(const std::complex<T>& z) {
	// 处理输入为零的情况
	if (z == std::complex<T>(0, 0)) {
		return std::complex<T>(0, 0);
	}

	// 计算复数的模
	T magnitude = std::abs(z);

	// 返回归一化的复数
	return std::complex<T>(z.real() / magnitude, z.imag() / magnitude);
}

vector<double> wavelet(vector<double> BLOCK1, int level)
{
	vector<double> swt_output;

	int length = 0;

	swt(BLOCK1, level, "haar", swt_output, length);

	return swt_output;
}

void printSignal(vector<double> &sig, int length = 0)
{
	if (length == 0) length = sig.size();

	for (int i = 0; i < length; i++)
	{
		cout << sig[i] << " ";
	}

	cout << endl;
}

double calculateMean(vector<vector<double>>& yo, int start, int end, int column)
{
	double sum = 0.0;
	int count = 0;

	for (int i = start; i <= end; i++)
	{
		sum += yo[i][column];
		count++;
	}

	return sum / count;
}

vector<vector<int>> general(const vector<vector<int>> &w, int m, int n, int a, int b)
{
	cout << "in general" << endl;

	vector<vector<int>> f(m, vector<int>(n, 0));

	for (int x = 0; x < m; x++)
	{
		for (int y = 0; y < n; y++)
		{
			int x1 = (x + 1) + a * (y + 1);
			if (x1 >= m || x1 < 0)
			{
				x1 = (x1 % m + m) % m;
				if (x1 == 0) x1 = m;
			}
			f[x1 - 1][y] = w[x][y];
		}
	}

	vector<vector<int>> tmp = f;

	for (int x = 0; x < m; x++)
	{
		for (int y = 0; y < n; y++)
		{
			int y1 = b * (x + 1) + (y + 1);
			if (y1 >= n || y1 < 0)
			{
				y1 = (y1 % n + n) % n;
				if (y1 == 0) y1 = n;
			}
			f[x][y1 - 1] = tmp[x][y];
		}
	}

	return f;
}

vector<vector<int>> matToVector(const cv::Mat &ww, int rows, int cols)
{
	vector<vector<int>> ww_matrix(rows, vector<int>(cols, 0));

	for (int i = 0; i < rows; i++)
	{
		for (int j = 0; j < cols; j++)
		{
			if (ww.at<uchar>(i, j) > 0)
			{
				ww_matrix[i][j] = 1;
			}
			else {
				ww_matrix[i][j] = 0;
			}
		}
	}

	return ww_matrix;
}

void printIntSignal(vector<int>& signal, int length = 0)
{
	if (length == 0) length = signal.size();

	for (int i = 0; i < length; i++)
	{
		cout << signal[i] << "\t";
	}
	cout << endl;
}

void printMat(vector<vector<double>> mat)
{
	for (std::vector<double> sub_mat : mat) {
		for (double value : sub_mat) {
			cout << value << "\t";
		}
		cout << endl;
	}
	cout << endl;
}

void printIntMat(std::vector<std::vector<int>> mat)
{
	for (std::vector<int> sub_mat : mat) {
		for (int value : sub_mat) {
			cout << value << "  ";
		}
		cout << endl;
	}
	cout << endl;
}

vector<double> readWav(drwav& wav, int& wavelength, int& channels, int& Fs) {
	// 使用 dr_mp3 打开 MP3 文件
	//drwav wav;
	//if (!drwav_init_file(&wav, filename, nullptr)) {
	//	std::cerr << "Error opening WAV file." << std::endl;
	//	return {};
	//}

	// 获取音频信息
	channels = static_cast<int>(wav.channels);
	drwav_uint64 totalPCMFrameCount = wav.totalPCMFrameCount; // 获取总 PCM 帧数

	// 计算总样本数
	wavelength = static_cast<int>(totalPCMFrameCount);
	Fs = static_cast<int>(wav.sampleRate);

	// 创建存储样本的向量
	std::vector<float> audio_f(totalPCMFrameCount * channels);
	std::vector<double> audio(totalPCMFrameCount * channels);

	// 读取 PCM 数据
	drwav_read_pcm_frames_f32(&wav, totalPCMFrameCount, audio_f.data());

	// 关闭 WAV 文件
	drwav_uninit(&wav);

	// 将 float 转换为 double
	for (size_t i = 0; i < totalPCMFrameCount * channels; i++) {
		audio[i] = static_cast<double>(audio_f[i]);
	}

	return audio;
}

bool canFormSquareMatrix(const int wm_size) {
	if (wm_size < 0) return false; // 负数不能组成方阵
	int root = static_cast<int>(std::sqrt(wm_size));
	return (root * root == wm_size); // 检查是否为完全平方数
}

void getAudio20sFront(bool split, int max_wavelength, int Fs, int channels, int &wavelength, vector<double> yo_origin, vector<double> &yo)
{
	if (split)
	{
		yo.reserve(max_wavelength); // 预留空间以提高性能
		yo.assign(yo_origin.begin(), yo_origin.begin() + max_wavelength);
		wavelength = max_wavelength / channels;
	}
	else {
		yo = yo_origin;
	}
}

void splitAudioByTime(int max_wavelength, int Fs, int channels, int& wavelength, vector<double> yo_origin, vector<double>& yo, const double start_time, const double end_time, int start)
{
	yo.resize(max_wavelength);
	for (int i = 0; i < max_wavelength; ++i) {
		yo[i] = yo_origin[start + i];
	}
	wavelength = max_wavelength / channels;
}