#include <stdio.
h>
#include <stdlib.h>
#include <math.h>
#include <complex.h>
#define NUM_FREQ_BINS 256 // Number of frequency bins
#define NUM_FRAMES 100 // Number of frames in STFT
#define FRAME_SIZE 512 // Frame size (samples per frame)
#define HOP_SIZE 256 // Hop size (overlap)
#define ALPHA 0.98 // Smoothing factor
#define GAMMA_N 0.5 // Diffuse noise coherence
// Function prototypes
void compute_covariance(double complex X_left[NUM_FREQ_BINS][NUM_FRAMES],
double complex X_right[NUM_FREQ_BINS][NUM_FRAMES],
double complex cov_matrix[2][2], int k, int l);
double complex compute_noise_psd(double complex cov_matrix[2][2]);
double compute_prior_snr(double complex noise_psd, double complex X);
double compute_gain(double posterior_snr, double prior_snr);
void inverse_stft(double complex X_enhanced[NUM_FREQ_BINS][NUM_FRAMES], double
*output_signal);
int main() {
// Initialize matrices for STFT of left and right channels
double complex X_left[NUM_FREQ_BINS][NUM_FRAMES];
double complex X_right[NUM_FREQ_BINS][NUM_FRAMES];
// Initialize output matrix for enhanced STFT
double complex X_enhanced[NUM_FREQ_BINS][NUM_FRAMES];
// Array for storing time-domain reconstructed signal
double output_signal[NUM_FRAMES * HOP_SIZE + FRAME_SIZE] = {0};
// Variables for noise PSD and SNR
double complex noise_psd[NUM_FREQ_BINS][NUM_FRAMES];
double prior_snr[NUM_FREQ_BINS][NUM_FRAMES] = {0};
double posterior_snr[NUM_FREQ_BINS][NUM_FRAMES];
double gain[NUM_FREQ_BINS][NUM_FRAMES];
// Fill X_left and X_right with the STFT data (replace with actual data
reading)
// ...
// Process each frequency bin and frame
for (int k = 0; k < NUM_FREQ_BINS; k++) {
for (int l = 0; l < NUM_FRAMES; l++) {
// Step 1: Calculate covariance matrix
double complex cov_matrix[2][2];
compute_covariance(X_left, X_right, cov_matrix, k, l);
// Step 2: Compute Noise PSD from second eigenvalue
noise_psd[k][l] = compute_noise_psd(cov_matrix);
// Step 3: Compute SNR
posterior_snr[k][l] = cabs(X_left[k][l]) * cabs(X_left[k][l]) /
creal(noise_psd[k][l]);
prior_snr[k][l] = ALPHA * prior_snr[k][l] + (1 - ALPHA) * fmax(0,
posterior_snr[k][l] - 1);
// Step 4: Calculate gain using GMMSE
gain[k][l] = compute_gain(posterior_snr[k][l], prior_snr[k][l]);
// Step 5: Enhance signal by applying gain
X_enhanced[k][l] = gain[k][l] * X_left[k][l];
}
}
// Step 6: Inverse STFT and overlap-add to reconstruct time-domain signal
inverse_stft(X_enhanced, output_signal);
printf("Denoising and reconstruction complete.\n");
// Output the enhanced signal (e.g., save to file or process further)
// ...
return 0;
}
void compute_covariance(double complex X_left[NUM_FREQ_BINS][NUM_FRAMES],
double complex X_right[NUM_FREQ_BINS][NUM_FRAMES],
double complex cov_matrix[2][2], int k, int l) {
double complex phi_xx = X_left[k][l] * conj(X_left[k][l]);
double complex phi_yy = X_right[k][l] * conj(X_right[k][l]);
double complex phi_xy = X_left[k][l] * conj(X_right[k][l]);
cov_matrix[0][0] = phi_xx;
cov_matrix[0][1] = phi_xy;
cov_matrix[1][0] = phi_xy;
cov_matrix[1][1] = phi_yy;
}
double complex compute_noise_psd(double complex cov_matrix[2][2]) {
// Calculate eigenvalues of a 2x2 matrix
double a = creal(cov_matrix[0][0]);
double b = creal(cov_matrix[0][1]);
double d = creal(cov_matrix[1][1]);
double trace = a + d;
double det = a * d - b * b;
double lambda1 = (trace + sqrt(trace * trace - 4 * det)) / 2;
double lambda2 = (trace - sqrt(trace * trace - 4 * det)) / 2;
// Use the smaller eigenvalue as noise PSD
return lambda2 / (1 - GAMMA_N);
}
double compute_prior_snr(double complex noise_psd, double complex X) {
double power_x = pow(cabs(X), 2);
return fmax(0, power_x / creal(noise_psd) - 1);
}
double compute_gain(double posterior_snr, double prior_snr) {
// GMMSE gain
return prior_snr / (1 + prior_snr);
}
void inverse_stft(double complex X_enhanced[NUM_FREQ_BINS][NUM_FRAMES], double
*output_signal) {
for (int l = 0; l < NUM_FRAMES; l++) {
for (int k = 0; k < NUM_FREQ_BINS; k++) {
double magnitude = cabs(X_enhanced[k][l]);
double phase = carg(X_enhanced[k][l]);
double real = magnitude * cos(phase);
double imag = magnitude * sin(phase);
// Reconstruct time-domain signal using OLA
int start_idx = l * HOP_SIZE;
for (int n = 0; n < FRAME_SIZE; n++) {
output_signal[start_idx + n] += real * cos(2 * M_PI * k * n /
FRAME_SIZE)
- imag * sin(2 * M_PI * k * n /
FRAME_SIZE);
}
}
}
}
//////////////not sure if its working
#include <stdio.h>
#include <stdlib.h>
#define MAX_ROWS 1000 // Set the maximum number of rows based on your data
#define MAX_COLS 1000 // Set the maximum number of columns based on your data
int read_csv(const char *filename, double matrix[MAX_ROWS][MAX_COLS], int *rows,
int *cols) {
FILE *file = fopen(filename, "r");
if (!file) {
perror("Unable to open file");
return -1;
}
char line[4096]; // Buffer to store each line from the CSV
*rows = 0;
*cols = 0;
while (fgets(line, sizeof(line), file) && *rows < MAX_ROWS) {
int col = 0;
char *ptr = line;
while (*ptr) {
matrix[*rows][col++] = strtod(ptr, &ptr); // Convert to double and
move pointer
if (*ptr == ',') {
ptr++; // Skip comma
}
}
if (*cols < col) {
*cols = col; // Update column count if needed
}
(*rows)++; // Update row count
}
fclose(file);
return 0;
}
int main() {
double matrix[MAX_ROWS][MAX_COLS];
int rows, cols;
if (read_csv("data.csv", matrix, &rows, &cols) == 0) {
printf("Matrix loaded successfully. Rows: %d, Columns: %d\n", rows, cols);
// Print matrix (optional)
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
printf("%f ", matrix[i][j]);
}
printf("\n");
}
} else {
printf("Failed to load matrix from CSV.\n");
}
return 0;
}
/////////////
#include <stdio.h>
#include <stdlib.h>
#define MAX_ROWS 1000 // Adjust based on the size of your matrix
#define MAX_COLS 1000 // Adjust based on the size of your matrix
void write_csv(const char *filename, double matrix[MAX_ROWS][MAX_COLS], int rows,
int cols) {
FILE *file = fopen(filename, "w");
if (!file) {
perror("Unable to open file");
return;
}
for (int i = 0; i < rows; i++) {
for (int j = 0; j < cols; j++) {
fprintf(file, "%f", matrix[i][j]); // Write the element
if (j < cols - 1) {
fprintf(file, ","); // Add comma after each element except the
last one
}
}
fprintf(file, "\n"); // Newline at the end of each row
}
fclose(file);
printf("Matrix written to %s successfully.\n", filename);
}
int main() {
double matrix[MAX_ROWS][MAX_COLS] = {
{1.1, 2.2, 3.3},
{4.4, 5.5, 6.6},
{7.7, 8.8, 9.9}
};
int rows = 3;
int cols = 3;
write_csv("output.csv", matrix, rows, cols);
return 0;
}
