How to design a FIR filter by frequency response Datas for shaping white noise

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)')

Answer 1

[William Rose]

@venkat ta,

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.

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.

[venkat ta]

Thanks for the reply.. actually this choherence transfer function is based on measured voltage and current.. so I tried inverting the impedance spectrum (providing the analyser, the spectrum of TFs of two inputs and its wav files ) for signal design, but didn't work. so why the switch to coherence?

Your FIR filter works fine with very high orders, but crest factors are increased >1, so additional compressor tuning is required.

%% designing FIR filter coeffiecent from normalized datas
subplot(4,2,4)

fy1 = fy1(1:index_20kHz);
Cxy1 = Cxy1(1:index_20kHz);
fs = 40000;

A=1./Cxy1;
F=fy1/(fs/2);             % normalized frequency
semilogx(F,20*log10(abs(A)),'g','LineWidth',1);hold on

[index_1000Hz,jj]=find(abs(fy1-1000)<0.02);% finding 1kHz index 
A1= A(1:index_1000Hz-1);A2=A(index_1000Hz:end);
[A2,window] = smoothdata(A2,'gaussian'); A=cat(1,A1,A2);
semilogx(F,20*log10(abs(A)),'r','LineWidth',1);hold on

N = 2048;                % filter order
d = fdesign.arbmag(N,F,A);  % A is linear level
Hd = design(d);         % design FIR filter

%  fvtool(Hd)
[H, w] = freqz(Hd);

% clf
semilogx(w/pi,20*log10(abs(H)),'b-','LineWidth',1);set(gca,'FontSize',20);fontname(gca,"Courier");
xt=xlabel('$Normalized  Frequency (\times\pi rad/sample)$');set(xt,'Interpreter','latex');yt=ylabel('Magnitude (dB)');set(yt,'Interpreter','latex');
% title('Frequency response of digital filter')

grid on
grid minor
ylim([-0.25 max(20*log10(abs(H)))])
LEGEND={['Normalized Freq of ', '$COH^{^-1}$',' No smooth from 1kHz'],['Normalized Freq of ', '$COH^{^-1}$',' with smooth from 1kHz'],['FIR response from Normalized  ', '$COH^{^-1}$']};legA=legend(LEGEND,'location','NW','FontSize',15,'LineWidth',0.5);legA.NumColumns = 1;set(legA,'Interpreter','latex')
% box off

%% shaping the signal x by FIR filter
subplot(4,2,[5 6])

y = filter(Hd,x);

plot(t,x,'r','LineWidth',1)% original signal
hold on
plot(t,y,'b-','LineWidth',1) % signal shaped
xt=xlabel('Time (sec)');yt=ylabel('0 dBFS');set(yt,'Interpreter','latex');set(xt,'Interpreter','latex')
set(gca,'FontSize',20);fontname(gca,"Courier");
grid on
grid minor
maX=max(abs([(min([min(x) min(y)])) (max([max(x) max(y)]))]));
ylim([-maX maX])
LEGEND = {['original signal',' max: ' num2str(max(x)),' min: ' num2str(min(x))], ['Filtered (FIR)',' max: ' num2str(max(y)),' min: ' num2str(min(y))]};
legA=legend(LEGEND,'location','SE','Orientation','horizontal','FontSize',15,'LineWidth',0.5);legA.NumColumns = 1;set(legA,'Interpreter','latex');legA.NumColumns = 2;

[William Rose]

@venkat ta,

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.