![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1279940/image.png)
How to design a FIR filter by frequency response Datas for shaping white noise
7 次查看(过去 30 天)
显示 更早的评论
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1279900/image.png)
Coherence data frequency curve, so I need to design the FIR to match that subplot 4's curve..
it seems filter: fdesign.arbmag has problem to process the data, can you clear me to go futher ?
clear all; close all; clc;
%% coherance data import for designing FIR filter (for shaping the white noise signal)
FileName = 'Coherence of white-noise_48k_2 with EQ No shaped.mat';
load(FileName)% if file is in local directory
%% process the loaded mat file
fy1 = Lin_Coh1(:,1);
Cxy1 = Lin_Coh1(:,2);
%% The signal need to be shaped in the end of the code by fIR filter
FileName_x = 'white-noise_48k_2 with EQ No shaped.mat';
load(FileName_x)
N = length(x);Fs= 48000;
t = (0:N-1)'/Fs; % t : time axis
%%
[index_20kHz,jj]=find(abs(fy1-20000)<0.02);% finding 20kHz index
subplot(3,2,1)
semilogx(fy1(1:index_20kHz),Cxy1(1:index_20kHz), 'LineWidth',2); hold on % plot the Coherence 0 - 1 (0 in dB means signal consitency good)
xlabel('Frequency (Hz)');ylabel('Coherence');set(gca,'FontSize',20);
grid on
grid minor
LEGEND={'coh : 0 = poor, 1 = good'};
legA=legend(LEGEND,'location','SW','FontSize',15,'LineWidth',0.5);
legA.NumColumns = 1;set(legA,'Interpreter','none')
clear LEGEND
%%
subplot(3,2,2)
semilogx(fy1(1:index_20kHz),mag2db(Cxy1(1:index_20kHz)), 'LineWidth',2); hold on % plot the Coherence 0 - 1 (0 in dB means signal consitency good)
semilogx(fy1(1:index_20kHz),mag2db(1./Cxy1(1:index_20kHz)), 'Color',"#EDB120",'LineWidth',2); % plot the inverse Coherence to process the filter shaping : EQ profiling)
xlabel('Frequency (Hz)');ylabel('Coherence in dB');set(gca,'FontSize',20);
grid on
grid minor
LEGEND={'coh in dB' '1/coh in dB'};
legA=legend(LEGEND,'location','SW','FontSize',15,'LineWidth',0.5);
legA.NumColumns = 1;set(legA,'Interpreter','none')
clear LEGEND
%% Inverse coh has to be processed so it's highlighted alone in subplot
subplot(3,2,3)
semilogx(fy1(1:index_20kHz),mag2db(1./Cxy1(1:index_20kHz)),'Color',"#EDB120", 'LineWidth',2); % plot the inverse Coherence to process the filter shaping : EQ profiling)
xlabel('Frequency (Hz)');ylabel('Inv coherence in dB');set(gca,'FontSize',20);
grid on
grid minor
LEGEND={'Inverse coh has to be processed'};
legA=legend(LEGEND,'location','NW','FontSize',15,'LineWidth',0.5);
legA.NumColumns = 1;set(legA,'Interpreter','none')
%% Inverse coh x and y axis normalized for designing FIR filter
subplot(3,2,4)
Nor_freq = 2.*fy1(1:index_20kHz)./48000;
Nor_gain = (1./Cxy1(1:index_20kHz)) ./ (max(1./Cxy1(1:index_20kHz)));
semilogx(Nor_freq,mag2db(Nor_gain),'Color',"#EDB120",'LineWidth',2); % plot the normalized - inverse Coherence to process the filter shaping : EQ profiling)
xlabel('Norm Frequency (x pi rad/sample)');ylabel('Norm Coh in dB');set(gca,'FontSize',20);
grid on
grid minor
LEGEND={'Normalized from subplot3'};
legA=legend(LEGEND,'location','NW','FontSize',15,'LineWidth',0.5);
legA.NumColumns = 1;set(legA,'Interpreter','none')
clear LEGEND
%% designing FIR filter coeffiecents from normalized datas
N = 150;
B = 3;
F = Nor_freq;
A = Nor_gain;
d = fdesign.arbmag('N,B,F,A',N,B,F,A);
Hd = design(d);
% fvtool(Hd)
[H, om] = freqz(Hd);
% Display frequency response (normalized frequency)
subplot(3,2,5)
clf
plot(Nor_freq, abs(H))
xlabel('Normalized frequency (cycles/sample)')
title('Frequency response')
box off
%% shaping the signal x by FIR filter
subplot(3,2,6)
y = filter(b, a, x);
clf
plot(t,x)% original signal
hold on
plot(t,y) % signal shaped
legend('original signal', 'Filtered (FIR)')
xlabel('Time (sec)')
0 个评论
采纳的回答
William Rose
2023-1-31
I think you want to create an FIR filter whose response is equal to the invese of the coherence. I will assume this is what you want. You did not state this clearly, so I could be mistaken.
fdesign.arbmag() is not happy with very long vectors of frequencies and amplitudes which you have passed. Therefore I recommend that you pick a few points from your plot of inverse coherence versus frequency, and pass those values to fdesign.arbmag.
From inspecting your plot, I estimated frequency and amplitude (in dB) at 17 points.
fs=44200; % sampling frequency (Hz)
f=[0,5,7,10,20,50,65,70,100,200,500,1000,2000,4000,6000,8000,10000,fs/2];
gaindB=[0,0.1,0.1,0.1,0.1,0.3,0.62,0.5,0.1,0.02,0.03,0.07,0.16,0.35,0.65,1.1,1.5,3.5];
A=db2mag(gaindB); % gain
F=f/(fs/2); % normalized frequency
N = 150; % filter order
d = fdesign.arbmag(N,F,A);
Hd = design(d); % design FIR filter
fvtool(Hd)
The script above generates a filter with approximately the deisred characteristics. See plot below.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1279940/image.png)
The dashed brown line is the desired response (F,A). The blue line is the actual filter response. Even with N=150, which is a relatively high filter order, the filter is not able to accurately reproduce the small bump in the amplitude plot that occurs at 65 Hz.
Why do you want a filter that is the inverse of the coherence plot? That seems odd to me. If you want to whiten signal, you would use the inverse of the signal amplitude, not the inverse of the coherence.
Good luck.
2 个评论
William Rose
2023-2-7
I see you have made a FIR filter, Hd, of order 2048. That is a very high order! I would not use such a high order filter because of fears about overfitting, but obviously you are comfortable with it. The filter's magnitude is A =1/Cxy, where Cxy = the coherence spectrum, smoothed above 1 kHz. Cxy was determined elsewhere, using data not shown here. Then you filtered a signal x with filter Hd to obtain signal y. I do not know if x is one of the signals used to compute Cxy, or of x is the magnitude of the impedance spectrum associated with Cxy.
I didn't know you could do subplot(4, 2, [5 6]) to use two panels at once in a subplot. Thanks for the tip!
Good luck with your signal processing.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Digital Filter Analysis 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!