need help with this

turn all the amplitudes of frequencies which are higher than 3500 Hz into 0’s and amplify the rest of the spectrum (by just multiplying the spectrum by some large number so that the amplitudes are around 2000-3000 in magnitude

how do i do this ? like what would the code be like?

Answer 1

[Sajeer Modavan]

a = [10;15;3600;4500;45;67;100];

[ind,~] = find(a>=3500);
b = a;
b(ind) = 0;
c = b*30

[dulanga]

hey i tried this but i am still getting alot of interferences

the file to filter

https://www.dropbox.com/s/xg77p5hbi5qr79f/unfiltered_sound.wav1?dl=0

[y,fs]=audioread('unfiltered_sound.wav1');
f=fft(y);
[ind,~] = find(f>=3500);
b = f;
b(ind) = 0;
c = b*30;
load handel.mat
filename = 'handel.wav';
audiowrite(filename,c,fs);
clear c fs

[Walter Roberson]

find(abs(f)>=3500)

And you have to ifft.

Loading handel at the end is not doing anything useful for you.

[dulanga]

yes i tried this before too doesnt work i still get high pitch noises

[y,fs]=audioread('unfiltered_sound.wav1');
f=fft(y);
[ind,~] = find(y>=3500);
b = y;
b(ind) = 0;
c = b*3000;
x=ifft(c);
load handel.mat
filename = 'handel.wav';
audiowrite(filename,c,fs);
clear c fs

[Walter Roberson]

Ah, the task is not about 3500 magnitude, it is about 3500 Hz. You need to figure out from fs which array locations in the result of fft() correspond to 3500 Hz or more.

When you do, be careful, because you also need to change the slots for the complex conjugates of the entries. For example, with a small number of samples, the order might look like

0 1 2 3 4 5 6 7 8 9 -9 -8 -7 -6 -5 -4 -3 -2 -1

where here -9 stands for the complex conjugate of what is in the slot marked 9.

Notice that the first entry in the array, the one marked 0, does not occur at the end. The first entry corresponds to the total of the input, which is equal to the mean multiplied by the number of samples; input signals that had a true mean of 0 would have 0 in that slot.

[dulanga]

how do u find that out? Is it like this?

[y,fs]=audioread('unfiltered_sound.wav1');
f=fft(y);
[ind,~] = find(fs>=3500);
b = fs;
b(ind) = 0;
c = b*3000;
x=ifft(c);
load handel.mat
filename = 'handel.wav';
audiowrite(filename,y,fs);
clear c fs

[Walter Roberson]

No.

L = size(y,1);
f = fs * (0:L/2)/L;

and you have to mirror f(2:end-1) or f(2:end) to get the second half of the frequencies. (Whether to use end-1 or end depends on whether you have an even number of slots or an odd number of slots.)

[dulanga]

so i should get the fft of f instead of fs ? and what do u mean by mirroring exactyl ?

[Walter Roberson]

[y,fs]=audioread('unfiltered_sound.wav1');
f=fft(y);
L = size(y,1);
f = fs * (0:L/2)/L;

Now here finish mirroring f before proceeding

[ind,~] = find(f>=3500);
b = f;
b(ind) = 0;
c = b*3000;
x=ifft(c);

filename = 'filtered_sound.wav';
audiowrite(filename,y,fs);

About mirrororing:

With a small number of samples, the slot order might look like

0 1 2 3 4 5 6 7 8 9 -9 -8 -7 -6 -5 -4 -3 -2 -1

where here -9 stands for the complex conjugate of what is in the slot marked 9.

In terms of frequency, if the original signal was at 1 Hz, the frequency list might look like

[0 Hz, 2 Hz, 4 Hz, 6 Hz, 8 Hz, 10 Hz, 12 Hz, 14 Hz, 16 Hz, 18 Hz, 18 Hz, 16 Hz, 14 Hz, 12 Hz, 10 Hz, 8 Hz, 6 Hz, 4 Hz, 2 Hz]

If you were given the vector

[0 2 4 6 7 10 12 14 16 18]

then you need to "reflect" it to get the second half -- the [18 16 14 12 10 7 6 4 2]

You should be able to work out how to do that.

But remember:

• you never "mirror" the first slot, which is for the DC component (0 Hz)
• After you attach the mirrored values, the total length of the frequency list has to be the same length as the fft
• the last entry in the half list might not get reflected. The other entries except the first always are. Consider the difference between [0 1 2 3 3 2 1] -- highest frequency gets duplicated when length of the signal is odd; [0 1 2 3 2 1] -- highest frequency does not get duplicated when length of the signal is even

[Sajeer Modavan]

values of 'f' doesn't go above 3500, so nothing replaced. Also, why are you multiplying with 30, I was showing an example, in your case it will be 3000/max(b), so that you will get maximum of 3000 in magnitude. I don't know much about noise removal, sorry if I am wrong.

[dulanga]

i am confused

f=fft(y);
L = size(y,1);
f = fs * (0:L/2)/L;

there are two fs now

[Walter Roberson]

Well, Yes, in that short sample data a that you posted, there are too few entries for a frequency as large as 3500. But in unfiltered_sound.wav1, there would be frequencies of 3500 or more if the fs was 7000 or more. In particular, if the original was sampled at 8000 Hz (telephone quality) then there would be frequencies present from 3500 to 4000.

[Walter Roberson]

[y, fs]=audioread('unfiltered_sound.wav1');
f=fft(y);
L = size(y,1);
Freqs = fs * (0:L/2)/L;

Now here finish mirroring Freqs before proceeding
Freqs = something for you to work out
    
%continuing
[ind,~] = find(Freqs>=3500);
b = f;
b(ind) = 0;
c = b*3000;
x=ifft(c);
filename = 'filtered_sound.wav';
audiowrite(filename,y,fs);

[Sajeer Modavan]

may be I am wrong, bcz as I said early, I am not expert in this domain, I didn't understood what you need between 2000-3000, magnitude or frequency, magnitude actually varying between -1500 to 1500 in this case, if you want change it, you can change d = c*3000/max(c);

[y,fs] = audioread('unfiltered_sound.wav1');
L = length(y);
Nf = ceil(L/2);
dxft = fft(y,Nf);
f = 0:fs/(2*Nf-1):fs/2;
figure
plot(f,dxft)

[~,ind] = find(f>=3500);
b = f;
b(ind) = [];
c = dxft;
c(ind) = [];

figure
plot(b,c)

[Walter Roberson]

Mohamed, you are correct that it would be best if the multiplier was determined from the signal. However, the instructions were,

"by just multiplying the spectrum by some large number so that the amplitudes are around 2000-3000 in magnitude"

so apparently the assignment expects the user to put in an arbitrary constant similar to your suggested 30, or what I seem to have slipped using as 3000, as long as the amplitudes come out in the right range. dulanga should apparently be plotting abs( c) to check the magnitudes and adjust the multiplier (the 30 or 3000) to give an acceptable range.

[dulanga]

walter is freqs like nqysuist or smthing ?

and i am confused with mirroring

the L is even so its is f(2:end-1) used but how

and you find in Freqs what is greater than 3500 but what happens to f where the fft is used ?

this was the guide given

[Walter Roberson]

discrete fourier transform theory says that you can take a function, f(x), and represent it as a sum of coefficients, each associated with a frequency -- essentially a sum of sine waves. The first slot of fft output is associated with "0 Hz", the DC constant. The second slot of fft output is for one cycle over the signal; the third slot is for two cycles over the signal; the third slot is for three cycles over the signal, and so on until you get to the middle slot, which for signal length L is associated with L/2 cycles per signal. The fact that you do not go on from there to (L/2+1), (L/2+2) and so on does have to do with Nyquist: Nyquist tells you that you cannot resolve anything more than L/2 cycles in a signal length L.

Once you get to the middle, then you start counting down again, possibly a repeat of L/2 (but not necessarily), and then L/2-1, L/2-2 and so on down to 1 (but not down to 0). When the input signal was entirely real-valued, then the coefficient in fft output array location K is the complex conjugate of the coefficient in location end-K+2. In the case of an input signal known to be entirely real-valued, you can reconstruct the entire signal from the first half of the fft() output, but it is easier to use the setup the way it is, as the current setup extends naturally to complex signals.

Now, when I say "one cycle over the signal, two cycles over the signal", and so on, I know that sounds like "Hz". But Hz is "cycles per second", and the signal is not necessarily one second long. If the signal were 2 seconds long, then "one cycle over the signal" would be 1 cycle in 2 seconds, which would be 1/2 Hz, for example. Two cycles over the signal would be 2 cycles in 2 seconds, 1 Hz. Three cycles over the signal would be 3 cycles in 2 seconds, 3/2 Hz. And so on. Notice that the adjacent slots are a constant frequency apart from each other. To know how many Hz each slot in the fft output represent, you need to do a small calculation. The calculation can be done if you know the total time that the signal represents, but it can also be done without that total time (directly) if you know the sampling frequency of the original signal, fs, which you get as the second output of audioread(). If you have more than one second of input, then the slots will correspond to frequencies less than fs apart; if you have less than one second of input, then the slots will correspond to frequencies more than fs apart.

For your task, you need to convert the bin numbers (cycles per signal) into absolute frequencies, because you are asked to filter based upon absolute frequencies. The Freqs calculation I posted earlier calculates the first half of that array, but does not continue on to count backwards down to lower frequencies again.

For mirroring: suppose you had an input vector x1 = [0 -1 2 -3 4] . Now, what expression would you write to get out [0 -1 2 -3 4 4 -3 2 -1] (case where original signal was odd length) ? What expression would you write to get out [0 -1 2 -3 4 -3 2 -1] (case where original signal was even length) ?