MATLAB Answers

Matlab - FFT/PSD Problem for Preemphasis

3 ビュー (過去 30 日間)
klaus ebert
klaus ebert 2015 年 4 月 19 日
編集済み: Youssef Khmou 2015 年 4 月 20 日
So my plan is the following (create "adaptive" pre/deemp. using Matlab 2015a):
  • 1. Read an Audiodata ([y,fs]) and generate white Noise with a certain SNR ([n,fs])
  • 2. Generate a Filter H which shapes the PSD(y) similiar to the PSD(n)
  • 3. Generate an inverse Filter G=H^(-1) which reverts the effect of H.
  • 4. Get the output signal
The code I was thinking to use for this is the following:
[y,fs]=audioread('test.wav');
snr=5;
n=awgn(y,snr,'measured')-y;
[pyy,f]=pwelch(y,[],[],length(y)*2-1,fs);
[pnn,fn]=pwelch(n,[],[],length(y)*2-1,fs);
H=sqrt(snr*pnn./pyy);
G=1./H;
outy=ifft(fft(y).*H.*G);
outn=ifft(fft(n).*G);
out=outy+outn;
[pout,fout]=pwelch(outy,[],[],length(outy)*2-1,fs);
[pnout,fnout]=pwelch(outn,[],[],length(outn)*2-1,fs);
There are two problems when using the code above:
  • 1. The pnout should be shaped liked pout (which it isn't).
  • 2. The outputsignals (outy,outn) are both complex and therefore it is not possible to play the outputfile.
If I try it another way:
[y,fs]=audioread('test.wav');
snr=5;
n=awgn(y,snr,'measured')-y;
N=length(y);
bin_vals=0:N-1;
fax_Hz= bin_vals*fs/N;
N_2=ceil(N/2);
pyy=(fft(y).*conj(fft(y)));
pnn=(fft(n).*conj(fft(n)));
H=sqrt(snr*pnn./pyy);
G=1./H
outy=ifft(fft(y).*H.*G);
outn=ifft(fft(n).*G);
out=outy+outn;
I don't come across the previous problems but in this case I don't use a real PSD but the squared absolute of the frequencydomain-values. (and as I understand it this would give me the power spectrum and not the power spectral density)
Any ideas why my first way is wrong or how to change my second way to get a real "proper" PSD and not just the PS?
Thanks! Klaus

  0 件のコメント

サインイン to comment.

回答 (2 件)

Youssef  Khmou
Youssef Khmou 2015 年 4 月 19 日
When using the reverse fft function, we take the real part, in one of the equations H.*G is not necessary since G=H.^-1, another method to generate noise is given in the following version, however it still needs scaling operation.
[y,fs]=audioread('test.wav');
snr=5;
sy=std(y);
sn=sy/snr;
n=sn*randn(size(y));
[pyy,f]=pwelch(y,[],[],length(y)*2-1,fs);
[pnn,fn]=pwelch(n,[],[],length(y)*2-1,fs);
H=sqrt(pnn./(pyy*snr));
G=1./H;
outy=(real(ifft(fft(y))));
outn=(real(ifft(fft(n).*G)));
out=outy+outn;
[pout,fout]=pwelch(outy,[],[],length(outy)*2-1,fs);
[pnout,fnout]=pwelch(outn,[],[],length(outn)*2-1,fs);
figure;
plot(fnout,pout,fnout,pnout,'r');

  1 件のコメント

klaus ebert
klaus ebert 2015 年 4 月 19 日
Thanks for your fast support! I thought about taking the real part as well but this gives me the peak of pnout at the wrong frequencies (e.g. with a monoton 1khz sound the peak of the outputnoise is at 2khz instead of 1khz)
Are there other probelms with my code or did I have a wrong thought at some point?

サインイン to comment.


Youssef  Khmou
Youssef Khmou 2015 年 4 月 19 日
I think the first question is not well interpreted, maybe it consists of generating noisy version of y instead of white noise, in this case, the result is sufficiently accurate :
[y,fs]=audioread('test.wav');
snr=5;
n=awgn(y,snr,'measured');
[pyy,f]=pwelch(y,[],[],length(y)*2-1,fs);
[pnn,fn]=pwelch(n,[],[],length(y)*2-1,fs);
H=sqrt(pnn./(pyy*snr));
G=1./H;
outy=real(ifft(fft(y).*H.*G));
outn=real(ifft(fft(n).*G));
out=outy+outn;
[pout,fout]=pwelch(outy,[],[],length(outy)*2-1,fs);
[pnout,fnout]=pwelch(outn,[],[],length(outn)*2-1,fs);
figure;
semilogy(fout,pout,fnout,pnout,'r');grid on;
xlabel('fr [Hz]');
ylabel('magnitude ');

  4 件のコメント

表示 1 件の古いコメント
Youssef  Khmou
Youssef Khmou 2015 年 4 月 19 日
from these tests, i think the problem is related to the method used to estimate the psd, especially the overlap issue and number of segments, let us try a last solution using standard psd function instead of pwelch, this requires a scaling of fft, try the following :
N=2000;
%[pyy,f]=pwelch(y,[],[],length(y)*2-1,fs);
[pyy,f]=psd(y,N);
%[pnn,fn]=pwelch(n,[],[],length(y)*2-1,fs);
[pnn,fn]=psd(n,N);
H=sqrt(pnn./(pyy*snr));
G=1./H;
outy=real(ifft(fft(y,N/2+1)));
outn=real(ifft(fft(n,N/2+1).*G));
out=outy+outn;
%[pout,fout]=pwelch(outy,[],[],length(outy)*2-1,fs);
[pout,fout]=psd(outy,N);
%[pnout,fnout]=pwelch(outn,[],[],length(outn)*2-1,fs);
[pnout,fnout]=psd(outn,N);
klaus ebert
klaus ebert 2015 年 4 月 20 日
Still the same problem as with pwelch. I am slowly really freaking out about the stretchfactor of 2 on the x-axis for the noise. Thank you for your intense support, do you maybe have any more ideas?
Youssef  Khmou
Youssef Khmou 2015 年 4 月 20 日
Here is another solution using two sided psd rather than one sided (psd or pwelch) in this case, we use the absolute value of fft :
% data
clear;
fs=1000;
ts=1/fs;
f=fs/4;
t=0:ts:0.1-ts;
y=sin(2*pi*f*t);
% altered part
snr=5;
sy=std(y);
sn=sy/snr;
n=sn*randn(size(y));
pyy=abs(fft(y));
pnn=abs(fft(n));
H=sqrt(pnn./(pyy*snr));
G=1./H;
outy=y;
outn=real(ifft(fft(n,length(y)).*G));
out=outy+outn;
[pout,fout]=psd(outy);
[pnout,fnout]=psd(outn);
figure;
semilogy(fout,pout,fnout,pnout,'r');grid on;
xlabel('fr [Hz]');
ylabel('magnitude ');

サインイン to comment.

サインイン してこの質問に回答します。


Translated by