How to filter all the frequencies from the audio signal.

Question

Justas Kvederis 2020-11-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/660578-how-to-filter-all-the-frequencies-from-the-audio-signal

评论： Mathieu NOE 2020-12-8

Hello guys, so I have this problem, I do have an audio signal, from two seperate sources with some noise. I need to seperate thoose into proper frequencies and print them out, how would I do that ? I already have the spectrum of the signal, that i acquired using FFT. But i get a lot of frequencies grouped in two spots, i don't get it how i get only the frequencies i need ? I added picture of the spectrum above.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Mathieu NOE 2020-11-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/660578-how-to-filter-all-the-frequencies-from-the-audio-signal#answer_554668

hello

2 gifts for you today (lucky man !)

a code for general spectral analysis
examples of digital filters you can use to eliminate unwanted frequency bands or keep (with a bandpass filter) frequencies of interest; once the filter is designed you can apply it on your signal using filter or filtfilt

hope it helps

36 个评论
显示 34更早的评论隐藏 34更早的评论

Mathieu NOE 2020-11-24

在 MATLAB Online 中打开

So this is an example with a notch filter centered on the first peak around 88 Hz (to remove)

of course you can apply whatever filter if your aim is different ... use this as an example

also there was an error when you do fft; you should only pick the first half of the fft points , the second half is only a mirror image (can be immediately recognized on your first picture)

look in documentation / fft / ('Single-Sided Amplitude Spectrum of S(t)')

so , the code :

clc
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1024*4;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% test data
Fs = 24000;
% t = 0:1/Fs:10-1/Fs;
fid = fopen('G4_5.dat','r')
[signal,samples] = fscanf(fid, '%12f');
fclose(fid);
 
%% notch filter section %%%%%%
% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
fc = 88;       % notch freq
wc = 2*pi*fc;
Q = 3;         % adjust Q factor for wider (low Q) / narrower (high Q) notch
                % at f = fc the filter has gain = 0
w0 = 2*pi*fc/Fs;
alpha = sin(w0)/(2*Q);
            b0 =   1;
            b1 =  -2*cos(w0);
            b2 =   1;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
            
% digital notch 
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
freq = linspace(fc-1,fc+1,200);
[gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
gd1db = 20*log10(gd1);
figure(1);
plot(freq,gd1db,'b');
title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
xlabel('Frequency (Hz)');
ylabel(' dB')
 
% now let's filter the signal 
signal_filtered = filtfilt(num1z,den1z,signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sensor_spectrum, freq] = pwelch(signal,w,NOVERLAP,NFFT,Fs);
sensor_spectrum_dB = 20*log10(sensor_spectrum);% convert to dB scale (ref = 1)
[sensor_spectrum_filtered, freq] = pwelch(signal_filtered,w,NOVERLAP,NFFT,Fs);
sensor_spectrum_filtered_dB = 20*log10(sensor_spectrum_filtered);% convert to dB scale (ref = 1)
figure(2),semilogx(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before notch filter','after notch filter');
xlabel('Frequency (Hz)');ylabel(' dB')

if you're happy with it please accept my answer

thanks

Mathieu NOE 2020-11-25

在 MATLAB Online 中打开

no problem

I could identify clearly 4 tones , but the fifth is where ? burried in the noise ?

from there we can use narrow bandpass filters (and put them in paralle) to keep the tones and reduce the noise

for the time being , this is just the code slightly modified with decimation to speed up things (as there is not much info in the signal above 1000 Hz)

clc
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1024;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% test data
Fs = 24000;
% t = 0:1/Fs:10-1/Fs;
fid = fopen('G4_5.dat','r')
[signal,samples] = fscanf(fid, '%12f');
fclose(fid);
% decimate
decim = 12;
signal = decimate(signal,decim);
Fs = Fs/decim;
 
%% notch filter section %%%%%%
% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
fc = 88;       % notch freq
wc = 2*pi*fc;
Q = 3;         % adjust Q factor for wider (low Q) / narrower (high Q) notch
                % at f = fc the filter has gain = 0
w0 = 2*pi*fc/Fs;
alpha = sin(w0)/(2*Q);
            b0 =   1;
            b1 =  -2*cos(w0);
            b2 =   1;
            a0 =   1 + alpha;
            a1 =  -2*cos(w0);
            a2 =   1 - alpha;
            
% digital notch 
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
freq = linspace(fc-1,fc+1,200);
[gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
gd1db = 20*log10(gd1);
figure(1);
plot(freq,gd1db,'b');
title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
xlabel('Frequency (Hz)');
ylabel(' dB')
 
% now let's filter the signal 
signal_filtered = filtfilt(num1z,den1z,signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sensor_spectrum, freq] = pwelch(signal,w,NOVERLAP,NFFT,Fs);
sensor_spectrum_dB = 20*log10(sensor_spectrum);% convert to dB scale (ref = 1)
[sensor_spectrum_filtered, freq] = pwelch(signal_filtered,w,NOVERLAP,NFFT,Fs);
sensor_spectrum_filtered_dB = 20*log10(sensor_spectrum_filtered);% convert to dB scale (ref = 1)
figure(2),plot(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before notch filter','after notch filter');
xlabel('Frequency (Hz)');ylabel(' dB')

Mathieu NOE 2020-11-25

在 MATLAB Online 中打开

code updated with multiple band pass filters implementation (butterworth)

if you find the missing frequency , you can add that to the list and run the code

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1024;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% test data
Fs = 24000;
% t = 0:1/Fs:10-1/Fs;
fid = fopen('G4_5.dat','r')
[signal,samples] = fscanf(fid, '%12f');
fclose(fid);
% decimate
decim = 12;
signal = decimate(signal,decim);
Fs = Fs/decim;
 
% %% notch filter section %%%%%%
% % H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
% 
% fc = 88;       % notch freq
% wc = 2*pi*fc;
% Q = 3;         % adjust Q factor for wider (low Q) / narrower (high Q) notch
%                 % at f = fc the filter has gain = 0
% 
% w0 = 2*pi*fc/Fs;
% alpha = sin(w0)/(2*Q);
% 
%             b0 =   1;
%             b1 =  -2*cos(w0);
%             b2 =   1;
%             a0 =   1 + alpha;
%             a1 =  -2*cos(w0);
%             a2 =   1 - alpha;
%             
% % digital notch 
% num1z=[b0 b1 b2];
% den1z=[a0 a1 a2];
% keep only signal content extracted by bandpass filters 
fc = [89 291 390 640];      % Hz
bandwith = 10;      % Hz
f_low = fc-0.5*bandwith;
f_high = fc+0.5*bandwith;
N = 4;
signal_filtered = zeros(size(signal));
for ci = 1:length(fc)
    [b,a] = butter(N,2/Fs*[f_low(ci) f_high(ci)]);
    freq = linspace(10,1000,100);
    [gd1,pd1] = dbode(b,a,1/Fs,2*pi*freq);
    gd1db = 20*log10(gd1);
    figure(1);
    plot(freq,gd1db);
    title(' Bandpass');
    xlabel('Frequency (Hz)');
    ylabel(' dB')
    hold on
    % now let's filter the signal and sum up with previous iteration
    signal_filtered = signal_filtered + filtfilt(b,a,signal);
end
hold off
% freq = linspace(fc-1,fc+1,200);
% [gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
% gd1db = 20*log10(gd1);
% 
% figure(1);
% plot(freq,gd1db,'b');
% title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
% xlabel('Frequency (Hz)');
% ylabel(' dB')
%  
% % now let's filter the signal 
% signal_filtered = filtfilt(num1z,den1z,signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before / after filters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sensor_spectrum, freq] = pwelch(signal,w,NOVERLAP,NFFT,Fs);
sensor_spectrum_dB = 20*log10(sensor_spectrum);% convert to dB scale (ref = 1)
[sensor_spectrum_filtered, freq] = pwelch(signal_filtered,w,NOVERLAP,NFFT,Fs);
sensor_spectrum_filtered_dB = 20*log10(sensor_spectrum_filtered);% convert to dB scale (ref = 1)
figure(2),plot(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before filter','after filter');
xlabel('Frequency (Hz)');ylabel(' dB')

Mathieu NOE 2020-12-1

在 MATLAB Online 中打开

ok

this is an attempt to make a better code - automatic also

but there will always be some situations where you have to adapt it to your signal

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% test data
Fs = 24000;
% t = 0:1/Fs:10-1/Fs;
fid = fopen('G4_5.dat','r')
[signal,samples] = fscanf(fid, '%12f');
fclose(fid);
% decimate
decim = 12;
signal = decimate(signal,decim);
Fs = Fs/decim;
samples = length(signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
n_max = fix(log2(samples));
NFFT = 2^(n_max-1);    % 2^n with n chosen so that NFFT is below length of signal (even after decimation) by factor > 2 so some averaging is doable
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);     % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sensor_spectrum, freq] = pwelch(signal,w,NOVERLAP,NFFT,Fs,'power');
sensor_spectrum_dB = 10*log10(sensor_spectrum);% convert to dB scale (ref = 1)
% sensor_spectrum_dB = sensor_spectrum_dB-min(sensor_spectrum_dB); % y axis shift so dB values are positive (min value = 0)
%% findpeaks
df = Fs/NFFT;
MAXW = 5*df;
MPD = MAXW;
MPP = 20;
[peaks,loc,W,P] = findpeaks(sensor_spectrum_dB,'MinPeakDistance',MPD,'MaxPeakWidth',MAXW,'MinPeakProminence',MPP);
freq_peaks = freq(loc);
figure(1),plot(freq,sensor_spectrum_dB,'b',freq_peaks,peaks,'*r');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(' dB')
 
% %% notch filter section %%%%%%
% % H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
% 
% fc = 88;       % notch freq
% wc = 2*pi*fc;
% Q = 3;         % adjust Q factor for wider (low Q) / narrower (high Q) notch
%                 % at f = fc the filter has gain = 0
% 
% w0 = 2*pi*fc/Fs;
% alpha = sin(w0)/(2*Q);
% 
%             b0 =   1;
%             b1 =  -2*cos(w0);
%             b2 =   1;
%             a0 =   1 + alpha;
%             a1 =  -2*cos(w0);
%             a2 =   1 - alpha;
%             
% % digital notch 
% num1z=[b0 b1 b2];
% den1z=[a0 a1 a2];
% keep only signal content extracted by bandpass filters 
fc = freq_peaks;      % Hz
bandwith = 10;      % Hz
f_low = fc-0.5*bandwith;
f_high = fc+0.5*bandwith;
N = 4;
signal_filtered = zeros(size(signal));
for ci = 1:length(fc)
    [b,a] = butter(N,2/Fs*[f_low(ci) f_high(ci)]);
    freq = linspace(10,1000,100);
    [gd1,pd1] = dbode(b,a,1/Fs,2*pi*freq);
    gd1db = 20*log10(gd1);
    figure(2);
    plot(freq,gd1db);
    title(' Bandpass');
    xlabel('Frequency (Hz)');
    ylabel(' dB')
    hold on
    % now let's filter the signal and sum up with previous iteration
    signal_filtered = signal_filtered + filtfilt(b,a,signal);
end
hold off
% freq = linspace(fc-1,fc+1,200);
% [gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
% gd1db = 20*log10(gd1);
% 
% figure(1);
% plot(freq,gd1db,'b');
% title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
% xlabel('Frequency (Hz)');
% ylabel(' dB')
%  
% % now let's filter the signal 
% signal_filtered = filtfilt(num1z,den1z,signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [sensor_spectrum, freq] = pwelch(signal,w,NOVERLAP,NFFT,Fs,'power');
% sensor_spectrum_dB = 10*log10(sensor_spectrum);% convert to dB scale (ref = 1)
[sensor_spectrum_filtered, freq] = pwelch(signal_filtered,w,NOVERLAP,NFFT,Fs,'power');
sensor_spectrum_filtered_dB = 10*log10(sensor_spectrum_filtered);% convert to dB scale (ref = 1)
figure(3),plot(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before filter','after filter');
xlabel('Frequency (Hz)');ylabel(' dB')

Mathieu NOE 2020-12-3

data acquisition.pdf

hi

decimate is usefull when the frequency range of interest ist much lower than Fs/2 (which is also called the Nyquist frequency).

if you don't want to have decimation, simply put 1 as value and there will be no more decimation;

and I may not agree about resolution; decimation does not reduce resolution, it reduces the amount of fft computation (therefore the computation time).

NB : the frequency resolution is the same if nfft and Fs are modified by the same factor (= decimation factor).

for the parameters like MAXW, MPD, MPP I had to try experimentaly, because I didn't find a mathematical robust way to link those parameters to other fft parameters. sorry about that.

for NFFT, NOVERLAP , you have to get more familiar with fft averaging techniques ; it's fairly common use in fft spectral analysis ; it's probably to technique that is mostly used in all vibration & noise FFT analysers.

if you google fft averaging, you will get lot of information about signal processing and fft techniques. You can also read vendor's brochures.

attached a usefull reminder about a few topics related to signal processing

Justas Kvederis 2020-12-7

在 MATLAB Online 中打开

Hello again, look I am going to admit, your code is wayyyy too difficult for me. And even if it wasn't that difficult I can't just reuse your code, so I am asking for the simplest approach to this.

My task is to find the starting and ending points when they are distributed like this.

Basically I have 4 sine waves mixed with noise, and now I need to find start point and end point to each of them. But I am really struggling with your code no matter how many times I look at it, for example I find my frequencies with simple for loop:

for i=1:length(Y)
    
    if Y(i) > 0        
        daugiau_uz = daugiau_uz + 1;  
        nuliai = 0;
    else
        nuliai = nuliai + 1;
    end
    
    if daugiau_uz > maziausias_plotis
        j = 0;
        sk = daugiau_uz;
        while sk >= 1             
            if Y(i - j) > 0
                klonas_blogu(i - j) = Y(i - j);
                sk = sk - 1;
            end
            j = j + 1;
        end    
    end
    if nuliai > 5
        daugiau_uz = 0;
    end
end    

This is only a snippet, but what I did i just made all spectrum values below 1000 = 0 (Y(Y<1000)=0;) and in loop I search for 5 zeros in a row in the spectrum, if I Find it I know it's either noise or sine wave, but since I know that sine wave will allways be thinner than noise, I just take 4 thinnest values and thoose are my frequencies. Well that's not the point I need something as simple as this, a loop a few simple functions, and that's it, right now I am attempting something like this:

piece = 500;
tro = zeros(1,N2);
for i=1:piece:length(signalas)-piece
    tro(1:piece) = signalas(i:i+piece-1);  
    tem = abs(fft(tro,piece));
    if tem(318) > 5
    
        disp(i);
T=(1/Fd)*piece;
delf=1/T;
daznAs=0:delf:Fd-delf;    
        stem(daznAs,tem, 'k.-');
    end
end

I know that this is incorrect, but what I tried to achieve is I split my signal piece by piece(size is 500) then i calculate spectrum for every piece and try to see if my desired frequency (this time 318) shows up, if it does I calculate the time when it show's up. How do I make it work, or do you know better solution ?

Mathieu NOE 2020-12-8

ok

hope your professor will be happy !

Mathieu NOE 2020-12-8

ok - not everything goes smooth and shinny as we would like too

I just hope you are not too much annoyed with signal processing !!

请先登录，再进行评论。

How to filter all the frequencies from the audio signal.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

36 个评论
显示 34更早的评论隐藏 34更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

How to filter all the frequencies from the audio signal.

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

36 个评论 显示 34更早的评论隐藏 34更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

36 个评论
显示 34更早的评论隐藏 34更早的评论