Techniques to reject noisy neural network input Suppose an artificial neural network is used to approximate a sine wave (shown in red in the graph below), given the linear input variable $x$ (scaled such that the ANN input is $x_{\rm nn}\in[-1;1]$).

The mean squared error between the target function (sine wave) and the trained ANN output is as shown approximately $4.8\times10^{-3}$, if the ANN structure is 1-8-1 (one input $x_{\rm nn}$, 8 hidden nodes and 1 output node).

Now suppose the ANN structure is changed to 2-8-1, where the second input is normally distributed random noise. The output is then as shown below. The MSE has increased to approximately $1.2\times10^{-2}$ (using the exact same noise sequence as during training) or slightly larger if a new noise sequence is used when testing the network.

The results above were obtained using standard gradient-descent type training with regularisation. All nodes in the hidden layer have $\rm tanh()$ activation functions and the output node is linear, as is standard with function approximation.

What methods are available to automatically remove the second input?
That is, training methods allowing the ANN to automatically discover the fact that the second input is noise, and therefore not a useful input. (Therefore excluding all methods where the human programmer must train several different ANN architectures with different input combinations, and evaluate all of them to decide which inputs to use and which to exclude.)

This question extends to larger ANNs where it is difficult to see beforehand which inputs are useful. Ideally, all weights associated with nonsense inputs should automatically be forced to zero.

It seems like your model is overfitting because you are not providing it with enough training data.

If I understood correctly, you want the neural network to learn $\hat{y} = \hat{f}(x,n)$ where $n \sim \mathcal{N}(\mu,\sigma)$ is noise and $x$ is in the interval $[-1,1]$ by minimizing
$$
\frac{1}{N}\sum_{i=1}^N |y_i - \hat{y_i}|^2 = \frac{1}{N}\sum_{i=1}^N |f(x_i,n_i) - \hat{f}(x_i,n_i)|^2
$$

And the function you are trying to learn is $f(x,n) = sin(\pi x)$, right?.

The network should learn to ignore the noise $n$ given that you provide it enough examples for training. That is, for different values of $x$, your training data should have a lot of examples with different values for $n$, all having the same target $sin(\pi x)$. If you show the network enough examples where the $n$ value is not useful for making the prediction, it should learn to minimize its effect on the output. The more training data you have, the better the network can learn to do this.


sampling frequency of GPS Now I am trying to design FIR based adaptive filter for rejection of Jamming in GPS device. (just self-learning purpose, Jamming is simple tone)

I've designed NOT-BAD performance adaptive filter for low sampling frequency. (~44100 Hz)

But when I use this adaptive filter for high sampling frequency (16 MHz), transient width of filter is greatly increased.
(I think it is sure because transient width is proportional to 1/degree)

Since it is self learning purpose, I just want to know how big sampling frequency of normal GPS signal for digital filtering is.

If it is not so high as above, I'll quit this design and concentrate on other issues(quantization error etc).

If sampling frequency of normal GPS signal is high as above, I'll study more about that issue.

So, my question is :

How high sampling frequency of normal GPS is ?

Thank you.

C/A code (the main civilian channel as of now) has a main lobe of around 2MHz and depending on how much processing you want to do, it is normally sampled in the low MHz range, which you are in. (you could pull it off with 2-5MHz).

Do know that the signal must be mixed down to a lower frequency before you sample it. The original signal is at 1.57GHz.

The frequency response of a BPSK signal is a sinc and therefore it's infinite. Most of the energy though is concentrated in the first null-to-null lobe.
Here in this picture you can see the a rectangle representing a BPSK pulse and its transform. The rectangle in this example lasts 200 samples and the sampling frequency is 1000. So the amplitude of the first lobe is 5Hz at (complex) baseband and 10Hz at (real) IF.
If we sampled only at 5Hz at baseband (or 10Hz at IF), that would be like applying a window of 5Hz in the frequency domain, in red in the second picture. The effect in the time domain would be a smoothing of the original rectangle and a loss of signal power.
Because GPS signals are autocorrelated in the receiver, the result of the correlation is also affected by the small sample rate. Instead of being a triangle (see 4th picture), what we get is a smooth sinc. The result is a weaker signal (with a loss in the order of 1.8dB) but potentially also a loss of precision in the computed pseudorange.
For GPS L1 C/A, the BPSK modulation has a frequency (aka chipping rate) of 1023 chips/ms. Therefore the first null is at 1.023MHz and the bandwidth of the first lobe is 2.046MHz at IF. A sampling rate of exactly 2.046Mhz would be sufficient to capture the energy from the first lobe.
That particular number (an integer multiple of the chipping rate) though would be a poor choice for other reasons.


Polar form of the Fourier transform of $\sin(t)$ I came across the Fourier transform of sin(t).

It ends up being a purely imaginary (dirac delta) impulse pair. But when considering the frequency domain representation of a signal, we consider the graphs of (i) the magnitude of the signal w.r.t. frequency and (ii) the phase (i.e. argument or angle) of the signal w.r.t. frequency.

In the case of $\sin(t)$, the phase can be seen to be 90 degrees, treating the delta functions in manner analogous to real numbers.

But how would I find the magnitude? Could someone please provide a reference to a textbook?

As another example, what would be the magnitude and phase of

$$\operatorname{FourierTransform}(\cos(t)) \cdot (4 + 2j) \text{ ?}$$

I read two question here. One is about the graphical representation of the Fourier transform of a sinusoidal function, and the other is about the magnitude (and phase) of such a Fourier transform. For ordinary functions, these two would be equivalent, but for Fourier transforms containing Dirac delta impulses they are not. The reason is that mathematically you can't compute the magnitude of a Dirac delta impulse. It is not a function, so there are no function values to take the magnitude of.

A standard graphical representation of the Fourier transform of the function $\sin(\omega_0t)$ would usually just be two arrows at $\omega=\omega_0$ and $\omega=-\omega_0$, one pointing up, the other one pointing down, where the $y$-axis is labeled by $j$ (which makes a separate phase plot unnecessary). An example of such a graph is shown here.

Even though you can't compute the magnitude of a Dirac impulse, you could write down a representation of the Fourier transform of a sinusoid consisting of a sum of Dirac impulses with positive areas and a phase term. For the function $f(t)=\sin(\omega_0t)$ such a representation would be

$$F(\omega)=\pi[\delta(\omega-\omega_0)+\delta(\omega+\omega_0)]e^{-j\frac{\pi}{2}\text{sign}(\omega)}\tag{1}$$

(Note that I use the non-unitary definition of the Fourier transform, which is different from the unitary version used by Wolfram Alpha. This is why there's a factor of $\pi$ in $(1)$ instead of $\sqrt{\frac{\pi}{2}}$.)

A similar representation of your example $A\cdot\mathcal{F}\{\cos(\omega_0t)\}$ with complex $A=|A|e^{j\phi_A}$, where $\phi_A$ is the argument of the complex number $A$, is given by

$$A\cdot\mathcal{F}\{\cos(\omega_0t)\}=|A|\pi [\delta(\omega-\omega_0)+\delta(\omega+\omega_0)]e^{j\phi_A}\tag{2}$$

Regarding to the magnitude and phase representation of delta functions,

http://www.dspguide.com/ch11/1.htm
http://www.math.hcmus.edu.vn/~tatuana/Xu%20Ly%20Anh/Ebook%20Others/CS450%20Slide%20Khac/ft-properties.seminar.pdf

Now, for the second example, you may need to develop that function as:

$4\cdot\cos(t) + 2j\cdot\cos(t)$ and then, you will end up with some deltas as well. This is the answer for the $2j\cdot\cos(t)$ one: http://www.wolframalpha.com/input/?i=fourier+transform+2jcos(x)

Note that you don't really need to expand the second example, but if you want to find the polar representation, you are going to do it anyway.


Power spectrum of a pure sine problems I am writing a C++ program to use to analyse oscillations and create a power spectrum.
I am using FFTW3.

To test the system, I am using a perfect sine function.

But no matter what parameters I try, I run into the same problem. I find the correct frequency, but also many 'overtones'.

First I assumed the problem was my way of chopping up the signal or how I calculate the running mean. But the problem is there when I just have all data in a single big window, bypassing my suspect code, still creating the pattern in the pictures.

I then assumed it was an aliasing problem, but even 1024 datapoints in a period, and they don't go away. They just lose shape and appear randomly. I assumed it was too little periods in my window. But more periods per window and they just become less frequent and stronger.
Then I assumed it may be spectral leakage, so I assumed I may need to use a windowing function like a Hamming, but that creates an odd slant. Is that the trade-off?

Or do I just set everything that is small enough to zero? I was just expecting a single strong point peak and the remainder to be so near zero, the log is -6 or lower everywhere else uniformly.

If I can't figure out how to get the right parameters for the best signal-to-noise ratio, I am not sure how I will achieve that when I start to work with the actual data.

I assumed it would be simple to pick a window that contains 32 periods with 32 discrete points in each period, use Welch's method and be done with it.

I am not an engineer who took many many classes on this kind of stuff. So I just don't have the proper training. I am working on a project of 4 weeks. There is so much info out there, it is confusing, but I don't know where to start as I don't have the time to learn stuff properly.

It can't be caused by the limited numbers of decimals of the input sine values, right?

Since FFT treats the signal as if it is periodic you need either to apply a window function (for example hanning) on your signal or make it coherent. In the image3 you attached you may make it coherent by only using data for a number of cycles; use samples 0..31 or 0..63. If the signal is non-coherent it will be seen by the FFT as a concatenated signal with discontinuites that will cause overtones.

EDIT:
A coherent signal may be concatenated with itself without generating any discontinuites.

The blue line is the input data to FFT

And a non-coherent:

I think I found the cause of the issue. It was the rounding off to three decimal places of the numbers I sampled from the sine. When I write output with 20 decimals, the overtones are gone. With 8 accuracy, they are still there but lower in amplitude. They also get more variance.

A new issue I wonder about is why my noise floor is -8 dB. Looking at other examples of power spectra created by DFT, I often see -20 or -50 dB for the frequency bins that are 'empty'.

Also, not sure if I can ignore this whole problem or if I have to apply a windowing function when it comes to the real data. I don't want aliases of the the strongest signal hiding weaker signals. And I don't know if I should like what the Hamming does. Maybe I am just irrationally biased against the slant it gives for the pure sinusoid.
Either way, I don't think putting in data with 20 decimal places is the right way to go, because that will make the input datafiles a lot bigger.

So the true reason was a pattern in the rounding of the sine created. If I copy&paste the first period, the FFTW3 library actually returns exact zero's for all frequency bins except integer multiples of the principle frequency. Those multiples are -13 dB and the signal itself is 3dB.
Without manually copy&pasting, a few points of the sine I generate are off by value of magnitude 10^-13 of what they were the period previously.

But I can't get any scalloping/spectral leakage as is described what one should get from the sinc function. I try to have non-integer number of periods in a window. I try to have the periods off-center. But I don't get actual scalloping. Just an uniform noise floor. Only thing kind of resembling scalloping I got from bad rounding. But those were more like aliases, right. Not an artifact of the sinc function.


How to add awgn noise to a signal in Fourier domain for given SNR I want to know the equations required to add awgn noise to a signal in Fourier domain for different SNR values. I know the procedure in time domain but not able to do the same in  frequency domain for the signal $x(t) = e^{2 \pi j f t}$ .

Thank you

Since the Fourier transform (and its discrete version) is linear, there is literally no difference in the way you do this in time and frequency domain:

You observe the energy in your signal, set the Gaussian's variance according to the desired SNR, and add the noise. Not a single flop was given about in which domain you were that day!


FFT analysis of binned data of different length I am conducting experiments to collect wind speed data from wind anemometers placed on a moving platform. Closely, there is a fixed wind mast holding a wind vane. Prior to analysing the wind data from the moving platform, I binned the data with 1-minute average wind direction. Each bin contains a number of continuous data sets (Set 1, Set 2, Set 3 etc.) which may not necessarily have the same number of data points. For example, Set 1 has 20,000 data points, Set 2 has 6,000 data points, etc.

I am performing an FFT on each of the sets to detect any periodicities in the signal. I am plotting the real amplitude of the sinusoids present in the signal on the y-axis vs. the frequency of the sinusoid. The FFT of each data set is yielding different amplitude results, on the y-axis and I'm not sure whether this is a result of the different length of the arrays. I am interested in an average value of the y-axis of the spectrum. However, I don't think that it is suitable to just calculate the average due to the different array size. Can someone please suggest a method? I am not used to analysis of signals in this way so I'm a bit confused. Thanks a lot.

Let's say you have several data sets $x_i$ and $N_i$ is the length of the $i$th data set.

First you're right that the amplitude of the FFT output scales with its length. To normalize all data sets they have to be divided by $N_i$. (Note that depending on the implementation of the FFT there might be a scaling involved already, check the Doc or see this question)

To combine all data sets you should append $(\max_i N_i) - N_i$ zeros to every data set before taking the FFT. The frequency axes are then aligned correctly and you can calculate the average of all FFT outputs.

Btw, I would plot the absolute value of the FFT output, not the real part. And sometimes I find it useful to "smooth" the spectrum by averaging every $M$ neighbouring FFT bins. Last, I have assumed that the data of every set has been acquired at equidistant time steps, otherwise the FFT isn't applicable at all.


Apply band pass filter (BPF) before calculating the power of a signal (PSD)? From my understanding PSD tries to give a good estimation on the power each frequency attributes to the overall signal power.

If I am only interested in a frequency range $(\omega_0 \leq \omega \leq \omega_1)$, could I put the signal through a band pass filter for that same frequency range then calculate the signal power using this formula:

$$P_x = \frac{1}{N_1 - N_0 + 1}\sum_{n = N_0}^{n=N_1} \left|x(n)\right|^2  $$

Would this technique give me what I am looking for, (the power of the signal from $N_0$ to $N_1$ for the specific frequency range)?

How does this compare to PSD estimation using any of the popular techniques?

Pros, Cons?

Thank you!

You can switch to the frequency domain with regards to Parseval's theorem. I guess you can transform your problem via discrete-time fourier transformation.


input output read/write frequency I am trying to identify a system based on input/output response and thereby estimate a transfer function.
I generated a frequency sweep function in Mathematica which gives me the discrete values of a frequency sweep.

g[t_] := 0.799760*Sin[2 Pi (3 t/333.3 + 1) t];
Plot[g[t], {t, 0, 1000000000000000000000000000000}]
g[#] & /@ Range[0, 5, 0.001]

So now I have a bunch of data but at what rate should I input it into the system, and at what rate should I read it?

How should I apply Nyquist frequency rule here?

Best to base it on the highest frequency response you anticipate for the system and implement twice that frequency with the samples you input and the data you read.


From pressure to SPL : microphone calibration using sound level meter calibrator and sound level meter I have a microphone + signal conditioner with adjustable gain and filter. I capture amplified signal using oscilloscope. I'm using this signal to compute equivalent continuous sound level

With signal conditioner, I set needed gain and A-weighted filter.
From oscilloscope signal V I can get pressure p = V/ (gain x mic sensitivity)

We can suppose A-weighted filter filters out frequency out of audio range. So I can easily compute leq = 20 log10(rms(pa)/pref), pa is time domain pressure A filtered, pref is reference pressure = 20e-6 pa

I used a sound level calibrator to check if it works. I get a leq of ~94dB and 114dB when applying 94dB and 114dB level with calibrator.

I also have a sound level meter, I also checked it with sound level calibrator.
Now, in my room (some noise because of oscilloscope...) with A filter it measures ~46 dB(A). Now I put my microphone close to it, with microphone I get a leq of 55 dB. When changing sound level meter filter from A to C results are better : sound level meter measures ~56dB.
I changed microphone acquisition duration but leq results are same.

How can we explain such leq difference when sound level meter and
microphone measure dB(A)?
How can I calibrate my microphone with    sound level meter?

First, I would check the positioning of sound level meter and microphone. If you are putting them close to the oscilloscope a few centimeters of displacement may change the measured value dramatically. (6dB/doubling of distance)

Second, which microphone do you use? The link you provided leads to a preamplifier. The actual capsule can either be a free field or a pressure field one. A free field capsule should be facing the sound source directly, while a pressure field capsule should face it sideways.


Calculating SNR for ECG signals I'm reading through the literature for finding the best way to calculate SNR for a ECG signal. While there are many algorithms for improving SNR, I couldn't find any that details any standard procedure for evaluating the SNR itself. The following are my questions.

What is the SNR of a typical ECG signal from the MIT/BIH database? What is considered a good SNR for a typical ECG signal?
http://www.physionet.org/physiobank/database/mitdb/
What would qualify as Signal and Noise for ECG? The assumption that I worked with is to apply a bandpass filter with cut-offs fc=[0.5Hz,45Hz] on the recorded ECG signal to get the Clean ECG signal, then subtract the clean ECG signal from recorded signal to obtain noise and then use the matlab snr function to evaluate  SNR with the clean signal and noise as arguments. But, in some papers it is recommended to reject everything except the QRS complex to get the clean signal. Is there any standard/recommended procedure for calculating SNR for ECG?
http://en.wikipedia.org/wiki/QRS_complex
What would be the best way to calculate ECG SNR in real-time in an embedded platform?

I don't know about real time application normally what I am doing is to calculate $\mathrm{PRD}$ first and calculate $\mathrm{SNR}$ from $\mathrm{PRD}$

$$\mathrm{SNR}=-20\log_{10}(0.01\mathrm{PRD})$$

$\mathrm{PRD}$ -  Percentage Root Mean Square Deviation


Time variant and Time Invariant Systems I was trying to solve the below system is a Time variant or Time Invariant System:

$y(n)=nx^{2}(n)$

Now for a system to be time invariant the system property should be :

y(n,k) = T[x(n-k)]

OR

y(n,k) = y(n-k)

for all possible values of k.

Coming back to system equation:

$y(n,k) = T[x(n-k)] = nx^{2}(n-k)$ delayed by k sample

And

$y(n-k) = nx^{2}(n-k)$

Hence

$y(n,k) = y(n-k)$

The the system is time invariant, but the solution in book states that the system is time variant. I am not sure where I am going wrong or my understanding is not correct?

Here is a simple way of looking at it which avoids any algebra and uses minimal symbolism. Let us calculate the output $y$ for an input $x$ which is a unit pulse function. We have that

$y[-2] = (-2)(x[-2])^2 = 0$ since $x[-2] = 0$.
$y[-1] = (-1)(x[-1])^2 = 0$ since $x[-1] = 0$.
$y[0] = (0)(x[0])^2 = 0$ (regardless of the value of $x[0]$).
$y[1] = (1)(x[1])^2 = 0$ since $x[1] = 0$.
$y[2] = (2)(x[2])^2 = 0$ since $x[2] = 0$.

and so on.  The entire input is blocked and the output $y$ is always $0$.

What about the output $\hat{y}$ when the input $\hat{x}$ is just the previous $x$
delayed by one time unit? The only difference now is that $\hat{x}[1] = 1$ and
we have that

$\hat{y}[1] = (1)(\hat{x}[1])^2 = 1$ since $\hat{x}[1] = 1$.

It is, I hope, obvious that $\hat{x}[n] = 0$ for all integers $n$ except $n=1$.
So, since the input $\hat{x}$ is just a delayed version of input $x$, is
is true that the output $\hat{y}$ is just a delayed version of output $y$?
What, if anything, does this tell you about the system properties?

$Y[n]=nx^2[n]$

$X[n]=nx^2[n]$

$Y[n-n_0] \ne X[n-n_0]$ because $Y[n-n_0]=n-n_0(x^2[n-n_0])$, and $X[n-n_0]=nx^2[n-n_0]$

Hence it's a time variant system.


Peak detection approach What are the peak detection algorithms in existence? I have noisy data, and I would like to implement peak detection for this data. Data is in reverse, actually I am trying to determine the bottom.

Here is a snapshot of the data in Excel. I like to detect both bottoms. I thought about passing the data through a low-pass filter and subsequently do a moving average where I determine the peaks and within the moving average I do another search. I have zero DSP background; this is just a common sense approach. I would like to hear what the experts recommend.

Ktuncer, there are a number of methods you can use here. One method that I would recommend is to use a Discrete Wavelet Transform, (DWT), and in particular, look at the Daubechies Wavelet. I would pick, say, Daub-14 / Daub-Tetra.

Basically what you really need to do is 'trend' your signal, and then from there, do a min or max pick. That will get rid of your outliers. A daub-14/daub-tetra wavelet transform can help you do this, and this helps especially since you do not know the nature of your signal. (Using daub-14, you can accurately represent polynomial signals of degree 14/2 = 7, and it looks like you wont need more than that).

The computation of this wavelet transform essentially 'compresses' your energy into a few indicies. Those indicies represent weights on basis vectors. The rest of the weights will (ideally) be near zero. When you have noise in your signal (as you do), those weights that were normally zero have some weights now, but you can simply zero them out and 'denoise' your signal. Once that is complete, you can then do a simple max/min detect.

There are more details involved, you may send me an email if you like to discuss how to implement it. I have done similar work on this before.

EDIT: Here are some images illustrating Daub-Tetra Denoiser:

I think one typical peak detection algorithm is like this where ref is peak(bottom).

for i=1,N   {
       if i=1   {  ref＝data(i)　｝
       else ｛ if　data(i)＜＝ref　｛ref ＝data(i) ｝｝
}


What is the real difference between DSP and AI/data science? I am very much into DSP, but lately everybody is talking about "Data science" and "deep learning", and from what I understand the whole thing about data science is to take a huge row of data, and extract meaningful data only from it.

Isn't it what you do with DSP, when you have a signal with noise and you extract the data from a raw of sampling? Or for example when you auto-correlate to find matches between signals?
Do you use the same statistical tools in both cases?

The purpose of both DSP and machine learning is to transform the input signal / data set into more meaningful information. This could be,

removal of noise from an audio recording
location of faces in an image
classification of objects in an image
etc

The main difference as I see it, is that in DSP the transform is designed by the engineer. The engineer will choose a set of signal processing operations that give the desired output. That choice is guided by experience and validated from the results of experiments.

In machine learning the transform is learnt. Typically, this requires a set of training signals with known outputs, over which the system will optimise its transform. The engineer still has to choose a machine learning architecture that can replicate the transform that is needed.

Often the two will be combined. And DSP is kind of machine learning anyway - the machine learning bit is done by the human and the test set is the sum of their existence.

DSP code that i see and try to write usually has a lot fewer conditional execution constructs (if statements) than AI code does.


How can I classifying this License Plate images? Matlab , LPR , I want to perform automatic character recognition on licence plates. One of my preprocessing step is binarization. However I have four different class of images and I have to use a different automatic thresolding algorithm on each class to get good binarization results.

I would like to choose the right algorithm automaticaly, my idea is to first perform classification on the licence plates images. In the following picture, each row represent a class.

picture:

http://www.mathworks.com/matlabcentral/answers/uploaded_files/31813/class.jpg

I'm working with MATLAB, how can I do that?? what is the best way?
Any ideas are appreciated.

Here's a procedure you should try in a first attempt :

You have RGB images, so each images is a 3*RC matrix, where R is the number of rows and C the number of columns.
We will consider each pixel of an image as a realization of a random vector $\mathbf{x}$ of dimention 3. Hence, having an image, you have RC realization of $\mathbf{x}$. With the classification exemple you gave, each classe has 4 images, so we have 4*R*C realization of $\mathbf{x}$. We are abble to estimate some stastical moments, that we hope will be good descriptor, such as mean, covariance etc ... (mean can be estimated with empirical mean or median).
Let's just stick with the mean for now. Let us call $m_1,m_2,m_3,m_4$ the means of each class you gave.
Now we have a new image we want to classify, we estimate it's 3-dimensional mean $\mathbf{\mu}$, then using nearest neighbour, the class is given by :
$$
C=\arg_i\min\left(\lVert\mathbf{\mu}-m_i \rVert _2\right)
$$

This is very basic.
You may improve this scheme using additional features like variance (diagonal of the covariance matrix) or covariance matrix. You may estimate variance or covariance using empirical variance or empirical covariance estimator (realy just the classical formulae). However the euclidian distance won't be a good distance since you will comparing inhomogenous vector (the features will be composed of mean + covariance), you'll want to use a more advanced distance such as Bhattacharyya distance for the special case of gaussian distribution.

Take a look at this example that uses MSER regions to detect text, and the ocr function to recognize it.


How to downsample a matrix in columns ? MATLAB I am studying wavelets in image processing. I would like to learn how to downsample a matrix in columns using MATLAB. I have used the downsample(x,n) command to downsample the given matrix in rows. But how to do the same in columns?

x =

 1     2     3
 4     5     6
 7     8     9
10    11    12

y = downsample(x,2)

y =

 1     2     3
 7     8     9

But I would like to get the following result.

y=
     1       3
     4       6
     7       9
    10      12

How to do it in MATLAB? Please help.

Maybe Off-Topic here, but anyway !

Try it

y = downsample(x',2)'

Do not use the downsample function.

I used simple method in MATLAB:

x= [1 2 3;
    4 5 6;
    7 8 9;
    10 11 12];

x1=x(:,1:2:end)

Have a good time. ^^


Impulse response of a system I have the following difference equation:
$y[n]=y[n-2]-0.25y[n-4]+x[n]-x[n-2]$
I need to find the first 6 values of the impulse response of the system characterized by that equation. When I do that, the values I get are:
$[1,0,0,0,-0.25,0]$, starting from index zero. I have the answer, which states that the right response is $[1,0,0,0,-0.25,1]$. I really cannot find my mistake; any help will be appreciated.

Answers in textbooks are not necessarily right.. I have got the same answer as you ;)

If you are not sure that your method is right, you can always test digital systems on your PC (matlab, scilab, python, octave, or almost in any language )


why is moving towards higher frequencies so important in frequency spectrum? As we know, as frequencies becomes higher, bandwidth becomes higher.And, according to channel capacity theorem, channel capacity increases with higher bandwidth.

Also,  energy is directly proportional to frequency(E=hf). So, if frequency increases, signals possesses higher energy and can travel far.

Which one of them would be more relevant to prove our need of higher frequency? Why not the other one?

So, if frequency increases, signals possesses higher energy and can
  travel far.

Not really, it's more the other way round:

When transmitting signals with radiowaves in free-space, most telecommunications engineers use the Friis equation which gives the amount of received power $(P_r)$ according to emitted power $(P_t)$, distance between emitter and receiver ($R$) and wavelength ($\lambda$) :
$$
\frac{P_r}{P_t} = G_tG_r\left(\frac{\lambda}{4\pi R}\right)^2
$$
$G_t$ and $G_r$ being the transmitter and receiver's antenna gain.
We see with this equation that $\frac{P_r}{P_t}$ is proportionnal to $\frac{1}{f^2}$ (because $\lambda = \frac{c}{f}$) which means that the higher the frequency, the higher the loss (and ^2 higher !).
The equation $E = h\nu$ is more related to quantum physics and has no practical use in telecommunications. Maybe it means that higher frequency signals are "more energetic" that low frequency ones, I don't really know honestly, maybe someone else* can answer better.

So the need to increase frequencies is mostly due to :

Saturation of the currently used spectrum : if you read international and national regulations and standards, you'll see that being able to emit on the RF spectrum is very expensive, because it's actually saturated, every frequency is reserved for a certain application (e.g. GSM signals is my country are transmitted in 200 KHz wide canals around 900 Mhz-ish frequencies). Still in my country, a new mobile phone operator had to pay billions of € to the governement to access some mobile comm. bands. Increasing the frequencies will give everyone more "room" to transmit radiowaves.
Increasing data rates, because having a higher carrier frequency also allows to widen the communication channel
Miniaturization (higher frequencies = lower wavelengths = smaller antennas and systems)

*edit :

  your claim on signal energy and frequency relation, which is based on photo-electricity (E=hv) effect is improper. First of all, the physical nature of waves has nothing to do in developing the mathematical notion of waves and their spectral representation by periodic signals. And even though it's true that the magnitude-squarred |x(t)|2 signal power relation is basicly adopted from power calculations in physics, that power (and energy) is implied to be the energy carried out by the waveform of a classical physics kind.

from Fat32's comment above

I'd like to make a comment about the relevance of this question in relation to 5G standards, where the lack of available spectrum is making researchers consider more and more to use of millimetre waves. This is a simple article about it: http://www.radio-electronics.com/info/cellulartelecomms/5g-mobile-cellular/millimetre-wave-links.php and this one is a bit more complete: http://spectrum.ieee.org/telecom/wireless/millimeter-waves-may-be-the-future-of-5g-phones note that this is from 2013. At the moment there are lots of publications talking about this topic and it is creating even a bit of controversy :D


How to mix two audio files recorded with different settings to obtain a realistic result I am a software engineer, currently working on a signal-processing related project, which I am new to. For evaluating the algorithm(which is written in R), I need to mix (merge) two different audio signals synthetically. The question is how to mix these two signals and obtain a “realistic” result?

To be more specific, consider the following scenario:

I have recorded Signal_A, and downloaded signal_B from the internet (therefore two different microphones/settings are used for recording).
The sound pressure level of both signals are known( I measured the sound pressure level of Signal_A in 1m distance, and also know the sound pressure level of Signal_B in the same distance).
As far as I know, the value recorded/represented in each waveform of any signal is relative to that corresponding signal, and since two signals are recorded using different settings we can expect that one might be more amplified than the other. Therefore a simple mixing of the two won’t give us a realistic result.

Now the question is: is there a way to manipulate recorded signal_A and signal_B programatically, assuming the fact that the SPL values of the original signals at 1m distance is known, and get the merged signal as they were playing simultaneously and recorded at 1m distance? i.e. to calculate factor ‘a’ and ‘b’ in a way that the realistic merged signal is formulated by:

Mereged_Signal = a*Signal_A  + b*Signal_B

Thanks in advance!

I have recorded Signal_A, and downloaded signal_B from the internet (therefore two different microphones/settings are used for recording)

You might be able to achieve a "seamless" mix (or transition) but the fact that you have two different microphones and settings might add artifacts in the final product. You can tell how this sounds like if you have ever noticed the differences between microphones. Accounting for such differences might be a bit more challenging.

  The sound pressure level of both signals are known( I measured the sound pressure level of Signal_A in 1m distance, and also know the sound pressure level of Signal_B in the same distance).

It is good that you have measured the SPL but most commonly you are unlikely to know the exact relationship that takes you from SPL to converted values.

  As far as I know, the value recorded/represented in each waveform of any signal is relative to that corresponding signal, and since two signals are recorded using different settings we can expect that one might be more amplified than the other. Therefore a simple mixing of the two won’t give us a realistic result.

I am not entirely clear on what you mean by "a simple mixing" but the first approach would be to simply adjust the amplitudes of the two recordings and then "add" them together as per your MergedSignal equation.

  Now the question is: is there a way to manipulate recorded signal_A and signal_B programatically, assuming the fact that the SPL values of the original signals at 1m distance is known, and get the merged signal as they were playing simultaneously and recorded at 1m distance?

Yes. BUT!, you have to make sure that the rest of the signal pathway was the same as well. Now:

Sound waves hit the transducer. The transducer converts them to voltage with a (hopefully) linear relationship of Volts per Pascal of pressure. In addition, this relationship is frequency dependent. But let's say that it is constant across the spectrum and call it $\alpha$.
Behind the transducer there is an amplifier whose job is to take the electrical signal and bring it up to some level. Again, this is achieved with a (hopefully) linear relationship. Let's keep things simple and assume that the job of this amplifier is to take the signal of the transducer and adapt it to a standard line level of approximately 2 Volts (peak-to-peak). Let's call this $\beta$.
The final stage is the Analog to Digital Converter (ADC). It converts the 2volts peak to peak to some range with a given word length (for example 0..255, -127..128 and so on).

So, assuming that $\alpha$ (mic), $\beta$ (amp/signal conditioning) are the same, you can indeed scale the recordings by the ratio of the SPLs. That is:

$$MergedSignal = (1.0 - \frac{SPL_A}{SPL_B}) \times Signal_A + \frac{SPL_A}{SPL_B} \times Signal_B$$

In the more general case and when you don't know the SPL, you can simply "match" the average level of the converted values. So, same relationship, instead of $SPL_A, SPL_B$ you have the averages of the signals.

Hope this helps.

Some remark to the last reply:

In my opinion you can even do the signal mixing with different $\alpha$ and $\beta$ if you consider the difference between the measured SPL value and the computed SPL value.

Let $SPL_{Am}$ be the measured SPL value of the signal A, let $SPL_{Ac}$ be the computed SPL value of the signal A (add the A-weighted squared amplitude values in the FFT frequency spectrum). Do the same for signal B. Then:

$$MergedSignal =  \frac{{SPL_A}_m}{{SPL_A}_c} \cdot  Signal_A + \frac{{SPL_B}_m}{{SPL_B}_c} \cdot Signal_B$$

Be careful with the dB-values. In the merging formula you need the non-logarithmic values: $10^{\frac{dB}{20}}$.

The difference between my formula and the one of A_A:

A_A's formula makes $Signal_A$ and $Signal_B$ equally loud. My formula makes $Signal_A$ as loud as it was (before scaling with $\alpha$ and $\beta$ due to microphone and amplifier), and $Signal_B$ as loud as it was.


Frequency response of a filter - fft vs. freqz I have designed my filter using firhalfband(N,Fo,A), which gives me the filter coefficients. Then I try to see the frequency response of it. I use these two function: fft, and freqz, but they are giving me different result in my figure, why is that? I wonder how these two MATLAB functions operate when taking the frequency response of a signal. I am expecting them to be the same.

Thank you

If you have the filter coefficients, freqz is the way to go - the various optional input arguments and the expected outputs are documented in MATLAB which you can see using help freqz. Snippets given below:

" freqz Digital filter frequency response.
    [H,W] = freqz(B,A,N) returns the N-point complex frequency response
    vector H and the N-point frequency vector W in radians/sample of
    the filter:

[H,F] = freqz(B,A,N,Fs) and [H,F] = freqz(B,A,N,'whole',Fs) return
    frequency vector F (in Hz), where Fs is the sampling frequency (in Hz).  "

If you have filter's impulse response - which will be different from the filter coefficients - you can use it to find the frequency response of the filter using the FFT.


Downsampling/Decimation - When is well done? I am working on decimation of signal, and I want to know which is the best way to understand if the downsampling is well done or not. At the moment I am comparing the FFT of the source signal and the downsampled signal and I observed a downward shift of it (I think it is due to the lesser quantity of samples), I had also a look into the time behavior. Do you have other advices/suggestions?

EDIT

From the comments/answers:

one method to understand if the downsampling process is well done is by comparing the FFT of the decimated signal and the original signal. The original signal must be attenuated by the decimation factor as suggested in the answer from @user2934229. Otherwise the signals can be normalized by the length of the signal (not suggested in my case because the signal is too long);

Ignoring your original question, and referring to your comment:

  I am observing the FFT of the downsampled signal presents an attenuation. I don't understand why.

The output of FFT (or DFT in general) isn't normalized - it is relative to the signal length (the total amount of samples being transformed). When decimating a signal (say, by a factor of M), you end up with a signal with 1/M length the original signal's length.

So, because DFT/FFT results are relative to signal length, a decimated signal with a 1/M length would have DFT/FFT components "attenuated" by a factor of M as well, even if the original signal was properly downsampled (i.e. bandlimited before being decimated).

To properly compare between the DFT of different signals, you should divide the output of DFT by the length of the transformed signal, to normalize the transform.

A decent ,,downsampling'' (let's say for obvious audio purposes) would, as you observed, just scale the spectrum and pad it with zeroes at the upper end.

The point would be that this is hard to do. Just inserting zeros (let's say one between each pair of samples) will mirror the spectrum into the upper half of the new spectrum. You would first insert zeros and then filter with as steep a filter as you can get. This will introduce error.

You were talking about decimation, which means removing samples. That would upsample your signal. I suppose you were referring to inserting samples instead.


Locating reference audio in longer audio I'm trying to solve the following problem:

Let's say I have a reference sound and a recording that I know contains at least one occurrence of the reference sound. How can I determine the positions in the audio where the sound was played? It's possible that the longer recording could have noise, and I cannot assume perfect time alignment with the sound that's playing in the recording (meaning, if I knew the content of the sound played in the recording other than the reference sound, I cannot just subtract that from the recording itself).

The method I'm experimenting with is as follows:

Grab windows the size of the reference sound from the recorded sound
For each window, compute the N dominant frequencies (by finding peak magnitudes in the FFT results)
Assign points to the dominant frequencies of the reference sound (where the most dominant gets a score of N, and the least gets a score of 1)
Compare the window's dominant frequencies to these and score it

The next few parts are just experimental, and used to narrow down results.

Do some filtering on the matches (for instance, a match with a score lower than half the possible score gets filtered out)
Look at overlapping matches and compute the cross correlation with the reference sound (using every possible time shift, prefer the maximum result). Prefer the match which requires the smaller shift
Perform that shift on the match's endpoints

The problem is that I want to also detect if part of the reference sound was played, and step 6 would not allow that if the reference sound was played two times in quick succession where one was cut off.

It seems like there has to be a solution that's not as ad-hoc as this. I'm very new to the field of DSP, so there may be something simple that I don't know about.

Any ideas?

Thanks!

It seems that a matched filter is good enough for my use case.


How Do I Measure the Time Duration of a Finite Length Discrete Sequence? Assume I have a five-sample time-domain sequence (none of the five samples are zero valued) and the time period between each pair of samples is one second. Measured in seconds, what is the time duration of my five-sample sequence?

The following graphic should clear up why a useful answer is 5 seconds, consistent with conversions between the discrete-time and continuous-time domains. A good example of this are discrete time approximations to continuous time integration or differentiation where the inclusion of $\Delta T$ (sample period) is made. With all such approximations, the time duration attributed to each sample is one sample interval. For the Forward Euler approximation specifically, each sample covers the time duration up to and not including the next sample. Such as in the first graphic the time duration for a sample at $t=0$ is one sample interval of one second; up to and not including the next sample at $t=1$ when sampled at one sample per second, with a total time duration for the one sample as one second.
The bottom figure then shows how this would apply to five samples starting at $t=0$ resulting in a total time duration of 5 seconds.

The plot below shows the Forward Euler discrete-time approximation for the integration of a continuous time waveform $x(t)$, which would be given as follows in the formula below and demonstrated in the graphic, best showing how each sample carries the time duration of the sample period (and as David has pointed out in the comments, this equally represents the zero-order hold approach in the output of a typical Digital to Analog converter):
$$\int_0^{NT}x(t)dt \approx \sum_{n=0}^{N-1}x[n]T$$

The actual range of time used for each sample is specific to the analog-digital translation approach used, but would cover one sample duration in all cases, leading to a total duration of 5 seconds for the OP's case.

It's $N/f_s$. Because consider what happens if it's $(N-1)/f_s$ instead; take a sequence $x$ of 6 samples, sampling once per second. Going by latter, defining "duration" as that between samples, we get

Now what is half the "duration" of $x$? Well, take half the sequence and measure the same:

Bad. If on the other hand we define duration as intervals midpointed by samples:

Good to go. This definition is also better grounded physically; suppose we wish to integrate. Doing so, we'd draw Riemann shapes (e.g. rectangles) in the discrete case, but regardless of shape, we'd have $N$ bins, so if finding e.g. average velocity, the physical interval spanned by integration would be $N/f_s$.
To be clear, this is largely a matter of definition, as a case can be made the other way; above's meant to show ways in which $N/f_s$ is more sensible.


What is meant by a system's "impulse response" and "frequency response?" Can anyone state the difference between frequency response and impulse response in simple English?

The impulse response and frequency response are two attributes that are useful for characterizing linear time-invariant (LTI) systems. They provide two different ways of calculating what an LTI system's output will be for a given input signal. A continuous-time LTI system is usually illustrated like this:

In general, the system $H$ maps its input signal $x(t)$ to a corresponding output signal $y(t)$. There are many types of LTI systems that can have apply very different transformations to the signals that pass through them. But, they all share two key characteristics:

The system is linear, so it obeys the principle of superposition. Stated simply, if you linearly combine two signals and input them to the system, the output is the same linear combination of what the outputs would have been had the signals been passed through individually. That is, if $x_1(t)$ maps to an output of $y_1(t)$ and $x_2(t)$ maps to an output of $y_2(t)$, then for all values of $a_1$ and $a_2$,

$$
H\{a_1 x_1(t) + a_2 x_2(t)\} = a_1 y_1(t) + a_2 y_2(t)
$$

The system is time-invariant, so its characteristics do not change with time. If you add a delay to the input signal, then you simply add the same delay to the output. For an input signal $x(t)$ that maps to an output signal $y(t)$, then for all values of $\tau$,

$$
H\{x(t - \tau)\} = y(t - \tau)
$$

Discrete-time LTI systems have the same properties; the notation is different because of the discrete-versus-continuous difference, but they are a lot alike. These characteristics allow the operation of the system to be straightforwardly characterized using its impulse and frequency responses. They provide two perspectives on the system that can be used in different contexts.

Impulse Response:

The impulse that is referred to in the term impulse response is generally a short-duration time-domain signal. For continuous-time systems, this is the Dirac delta function $\delta(t)$, while for discrete-time systems, the Kronecker delta function $\delta[n]$ is typically used. A system's impulse response (often annotated as $h(t)$ for continuous-time systems or $h[n]$ for discrete-time systems) is defined as the output signal that results when an impulse is applied to the system input.

Why is this useful? It allows us to predict what the system's output will look like in the time domain. Remember the linearity and time-invariance properties mentioned above? If we can decompose the system's input signal into a sum of a bunch of components, then the output is equal to the sum of the system outputs for each of those components. What if we could decompose our input signal into a sum of scaled and time-shifted impulses? Then, the output would be equal to the sum of copies of the impulse response, scaled and time-shifted in the same way.

For discrete-time systems, this is possible, because you can write any signal $x[n]$ as a sum of scaled and time-shifted Kronecker delta functions:

$$
x[n] = \sum_{k=0}^{\infty} x[k] \delta[n - k]
$$

Each term in the sum is an impulse scaled by the value of $x[n]$ at that time instant. What would we get if we passed $x[n]$ through an LTI system to yield $y[n]$? Simple: each scaled and time-delayed impulse that we put in yields a scaled and time-delayed copy of the impulse response at the output. That is:

$$
y[n] = \sum_{k=0}^{\infty} x[k] h[n-k]
$$

where $h[n]$ is the system's impulse response. The above equation is the convolution theorem for discrete-time LTI systems. That is, for any signal $x[n]$ that is input to an LTI system, the system's output $y[n]$ is equal to the discrete convolution of the input signal and the system's impulse response.

For continuous-time systems, the above straightforward decomposition isn't possible in a strict mathematical sense (the Dirac delta has zero width and infinite height), but at an engineering level, it's an approximate, intuitive way of looking at the problem. A similar convolution theorem holds for these systems:

$$
y(t) = \int_{-\infty}^{\infty} x(\tau) h(t - \tau) d\tau
$$

where, again, $h(t)$ is the system's impulse response. There are a number of ways of deriving this relationship (I think you could make a similar argument as above by claiming that Dirac delta functions at all time shifts make up an orthogonal basis for the $L^2$ Hilbert space, noting that you can use the delta function's sifting property to project any function in $L^2$ onto that basis, therefore allowing you to express system outputs in terms of the outputs associated with the basis (i.e. time-shifted impulse responses), but I'm not a licensed mathematician, so I'll leave that aside). One method that relies only upon the aforementioned LTI system properties is shown here.

In summary: For both discrete- and continuous-time systems, the impulse response is useful because it allows us to calculate the output of these systems for any input signal; the output is simply the input signal convolved with the impulse response function.

Frequency response:

An LTI system's frequency response provides a similar function: it allows you to calculate the effect that a system will have on an input signal, except those effects are illustrated in the frequency domain. Recall the definition of the Fourier transform:

$$
X(f) = \int_{-\infty}^{\infty} x(t) e^{-j 2 \pi ft} dt
$$

More importantly for the sake of this illustration, look at its inverse:

$$
x(t) = \int_{-\infty}^{\infty} X(f) e^{j 2 \pi ft} df
$$

In essence, this relation tells us that any time-domain signal $x(t)$ can be broken up into a linear combination of many complex exponential functions at varying frequencies (there is an analogous relationship for discrete-time signals called the discrete-time Fourier transform; I only treat the continuous-time case below for simplicity). For a time-domain signal $x(t)$, the Fourier transform yields a corresponding function $X(f)$ that specifies, for each frequency $f$, the scaling factor to apply to the complex exponential at frequency $f$ in the aforementioned linear combination. These scaling factors are, in general, complex numbers. One way of looking at complex numbers is in amplitude/phase format, that is:

$$
X(f) = A(f) e^{j \phi(f)}
$$

Looking at it this way, then, $x(t)$ can be written as a linear combination of many complex exponential functions, each scaled in amplitude by the function $A(f)$ and shifted in phase by the function $\phi(f)$. This lines up well with the LTI system properties that we discussed previously; if we can decompose our input signal $x(t)$ into a linear combination of a bunch of complex exponential functions, then we can write the output of the system as the same linear combination of the system response to those complex exponential functions.

Here's where it gets better: exponential functions are the eigenfunctions of linear time-invariant systems. The idea is, similar to eigenvectors in linear algebra, if you put an exponential function into an LTI system, you get the same exponential function out, scaled by a (generally complex) value. This has the effect of changing the amplitude and phase of the exponential function that you put in.

This is immensely useful when combined with the Fourier-transform-based decomposition discussed above. As we said before, we can write any signal $x(t)$ as a linear combination of many complex exponential functions at varying frequencies. If we pass $x(t)$ into an LTI system, then (because those exponentials are eigenfunctions of the system), the output contains complex exponentials at the same frequencies, only scaled in amplitude and shifted in phase. These effects on the exponentials' amplitudes and phases, as a function of frequency, is the system's frequency response. That is, for an input signal with Fourier transform $X(f)$ passed into system $H$ to yield an output with a Fourier transform $Y(f)$,

$$
Y(f) = H(f) X(f) = A(f) e^{j \phi(f)} X(f)
$$

In summary: So, if we know a system's frequency response $H(f)$ and the Fourier transform of the signal that we put into it $X(f)$, then it is straightforward to calculate the Fourier transform of the system's output; it is merely the product of the frequency response and the input signal's transform. For each complex exponential frequency that is present in the spectrum $X(f)$, the system has the effect of scaling that exponential in amplitude by $A(f)$ and shifting the exponential in phase by $\phi(f)$ radians.

Bringing them together:

An LTI system's impulse response and frequency response are intimately related. The frequency response is simply the Fourier transform of the system's impulse response (to see why this relation holds, see the answers to this other question). So, for a continuous-time system:

$$
H(f) = \int_{-\infty}^{\infty} h(t) e^{-j 2 \pi ft} dt
$$

So, given either a system's impulse response or its frequency response, you can calculate the other. Either one is sufficient to fully characterize the behavior of the system; the impulse response is useful when operating in the time domain and the frequency response is useful when analyzing behavior in the frequency domain.

Shortly, we have two kind of basic responses: time responses and frequency responses. Time responses test how the system works with momentary disturbance while the frequency response test it with continuous disturbance. Time responses contain things such as step response, ramp response and impulse response. Frequency responses contain sinusoidal responses.
Aalto University has some course Mat-2.4129 material freely here, most relevant probably the Matlab files because most stuff in Finnish. If you are more interested, you could check the videos below for introduction videos. I found them helpful myself.
I have only very elementary knowledge about LTI problems so I will cover them below -- but there are surely much more different kinds of problems!
Responses with Linear time-invariant problems

With LTI (linear time-invariant) problems, the input and output must have the same form: sinusoidal input has a sinusoidal output and similarly step input result into step output. If you don't have LTI system -- let say you have feedback or your control/noise and input correlate -- then all above assertions may be wrong. With LTI, you will get two type of changes: phase shift and amplitude changes but the frequency stays the same. If you break some assumptions let say with non-correlation-assumption, then the input and output may have very different forms.
If you need to investigate whether a system is LTI or not, you could use tool such as Wiener-Hopf equation and correlation-analysis. Wiener-Hopf equation is used with noisy systems. It is essential to validate results and verify premises, otherwise easy to make mistakes with differente responses. More about determining the impulse response with noisy system here.

References

Wikipedia article about LTI here

Very good introduction videos about different responses here and here -- a few key points below.


Reconstruction of audio signal from Spectrogram I have a set of songs for which I extracted the magnitude spectrogram using a Hamming Window with 50% overlap. After extracting the spectrogram, I did some dimensionality reduction using Principal Components Analysis (PCA). After reducing it to lower dimensionality, I reconstructed the spectrograms from lower dimensions. So now, there would be some error between the original spectrogram and the reconstructed spectrogram. I would like to convert this spectrogram back to the audio signal and play it, so that I would be able to know when reconstructed from lower dimensions, how does the audio sound.

Is there any function available in say Matlab. to convert a magnitude Spectrogram to an audio signal ??

If the spectrogram was computed as the magnitude of short time fourrier transforms from overlapping windows, then the spectrogram contains implicitly some phase information.

The following iterations do the job :

$$x_{n+1} = \text{istft}(S\cdot\exp(i\cdot\text{angle}(\text{stft}(x_n))))$$

$S$ is the spectrogram, $\text{stft}$ is the forward-short time Fourier transform, $\text{isft}$ is the inverse-short time Fourier transform.

Use Griffin-Lim algorithm to invert the audio signal from spectrogram, if you are not worried about computation complexity.


Why is the Fourier transform so important? Everyone discusses the Fourier transform when discussing signal processing. Why is it so important to signal processing and what does it tell us about the signal?

Does it only apply to digital signal processing or does it apply to analog signals as well?

This is quite a broad question and it indeed is quite hard to pinpoint why exactly Fourier transforms are important in signal processing. The simplest, hand waving answer one can provide is that it is an extremely powerful mathematical tool that allows you to view your signals in a different domain, inside which several difficult problems become very simple to analyze.

Its ubiquity in nearly every field of engineering and physical sciences, all for different reasons, makes it all the more harder to narrow down a reason. I hope that looking at some of its properties which led to its widespread adoption along with some practical examples and a dash of history might help one to understand its importance.

History:

To understand the importance of the Fourier transform, it is important to step back a little and appreciate the power of the Fourier series put forth by Joseph Fourier. In a nut-shell, any periodic function $g(x)$ integrable on the domain $\mathcal{D}=[-\pi,\pi]$ can be written as an infinite sum of sines and cosines as

$$g(x)=\sum_{k=-\infty}^{\infty}\tau_k e^{\jmath k x}$$
$$\tau_k=\frac{1}{2\pi}\int_{\mathcal{D}}g(x)e^{-\jmath k x}\ dx$$

where $e^{\imath\theta}=\cos(\theta)+\jmath\sin(\theta)$. This idea that a function could be broken down into its constituent frequencies (i.e., into sines and cosines of all frequencies) was a powerful one and forms the backbone of the Fourier transform.

The Fourier transform:

The Fourier transform can be viewed as an extension of the above Fourier series to non-periodic functions. For completeness and for clarity, I'll define the Fourier transform here. If $x(t)$ is a continuous, integrable signal, then its Fourier transform, $X(f)$ is given by

$$X(f)=\int_{\mathbb{R}}x(t)e^{-\jmath 2\pi f t}\ dt,\quad \forall f\in\mathbb{R}$$

and the inverse transform is given by

$$x(t)=\int_{\mathbb{R}}X(f)e^{\jmath 2\pi f t}\ df,\quad \forall t\in\mathbb{R}$$

Importance in signal processing:

First and foremost, a Fourier transform of a signal tells you what frequencies are present in your signal and in what proportions.

  Example: Have you ever noticed that each of your phone's number buttons
  sounds different when you press during a call and that it sounds the same for every phone model? That's because they're each composed of two different sinusoids which can be used to uniquely identify the button. When you use your phone to punch in combinations to navigate a menu, the way that the other party knows what keys you pressed is by doing a Fourier transform of the input and looking at the frequencies present.

Apart from some very useful elementary properties which make the mathematics involved simple, some of the other reasons why it has such a widespread importance in signal processing are:

The magnitude square of the Fourier transform, $\vert X(f)\vert^2$ instantly tells us how much power the signal $x(t)$ has at a particular frequency $f$.
From Parseval's theorem (more generally Plancherel's theorem), we have
$$\int_\mathbb{R}\vert x(t)\vert^2\ dt = \int_\mathbb{R}\vert X(f)\vert^2\ df$$
which means that the total energy in a signal across all time is equal to the total energy in the transform across all frequencies. Thus, the transform is energy preserving.
Convolutions in the time domain are equivalent to multiplications in the frequency domain, i.e., given two signals $x(t)$ and $y(t)$, then if

$$z(t)=x(t)\star y(t)$$
where $\star$ denotes convolution, then the Fourier transform of $z(t)$ is merely

$$Z(f)=X(f)\cdot Y(f)$$

For discrete signals, with the development of efficient FFT algorithms, almost always, it is faster to implement a convolution operation in the frequency domain than in the time domain.
Similar to the convolution operation, cross-correlations are also easily implemented in the frequency domain as $Z(f)=X(f)^*Y(f)$, where $^*$ denotes complex conjugate.
By being able to split signals into their constituent frequencies, one can easily block out certain frequencies selectively by nullifying their contributions.

  Example: If you're a football (soccer) fan, you might've been
  annoyed at the constant drone of the vuvuzelas that pretty much
  drowned all the commentary during the 2010 world cup in South Africa.
  However, the vuvuzela has a constant pitch of ~235Hz which made it
  easy for broadcasters to implement a notch filter to cut-off the
  offending noise.[1]

A shifted (delayed) signal in the time domain manifests as a phase change in the frequency domain. While this falls under the elementary property category, this is a widely used property in practice, especially in imaging and tomography applications,

  Example: When a wave travels through a heterogenous medium, it
  slows down and speeds up according to changes in the speed of wave
  propagation in the medium. So by observing a change in phase from
  what's expected and what's measured, one can infer the excess time
  delay which in turn tells you how much the wave speed has changed in
  the medium. This is of course, a very simplified layman explanation, but
  forms the basis for tomography.

Derivatives of signals (nth derivatives too) can be easily calculated(see 106) using Fourier transforms.

Digital signal processing (DSP) vs. Analog signal processing (ASP)

The theory of Fourier transforms is applicable irrespective of whether the signal is continuous or discrete, as long as it is "nice" and absolutely integrable. So yes, ASP uses Fourier transforms as long as the signals satisfy this criterion. However, it is perhaps more common to talk about Laplace transforms, which is a generalized Fourier transform, in ASP. The Laplace transform is defined as

$$X(s)=\int_{0}^{\infty}x(t)e^{-st}\ dt,\quad \forall s\in\mathbb{C}$$

The advantage is that one is not necessarily confined to "nice signals" as in the Fourier transform, but the transform is valid only within a certain region of convergence. It is widely used in studying/analyzing/designing LC/RC/LCR circuits, which in turn are used in radios/electric guitars, wah-wah pedals, etc.

This is pretty much all I could think of right now, but do note that no amount of writing/explanation can fully capture the true importance of Fourier transforms in signal processing and in science/engineering

the physics relevance of fourier transform is that it tells the relative amplitude of frequencies present in the signal . it can be defined for both discrete time and continuous  time signal.
Any signal can be represented as mixture of many harmonic frequencies. Fourier transform help in filter applications , where we need only certain range of frequencies then we first need to know what are the amplitudes of frequencies contains in the signal.


Negative dB signal to noise ratio Is it common to find systems having signal to noise ratio (SNR) less than 1 i.e. noise is more than signal? I am currently working with SNR of -5db to -20db. And how useful is it to work with sensor readings having negative SNR? I am asking it because the algorithm I am using to solve a problem works well with negative SNR but at positive and higher SNR, it is numerically unstable. This instability arises because of computation of hypergeometric function $_pF_q\left(;;\frac{z}{σ^2}\right)$. At high SNR, i.e. low $σ^2$, the argument $\frac{z}{σ^2}$ becomes very high, making the convergence of hypergeometric function very slow and the algorithm becomes numerically unstable.

Yes, GPS signals are below the noise floor. Link

It's quite common to work with signals that are below the noise. If you use signal of some length more than 1 you can exploit it to overcome noise impact so it's all right. As I've understood the issue is with algorithm convergence. It's unknown and it's difficult to define what the problem is (write some details if you're allowed). But maybe you should use $\frac{z}{\sigma^2+1}$ instead? So if SNR is pretty low, $\sigma^2$ is the major one and your algorithm will work well as you've written. And while SNR is high the argument will not be very high no more.

If you describe your approach briefly, somebody can probably tell you more about how to improve your algorithm.


Why look at power spectral density for stochastic processes? I have been told that for deterministic signals, it makes sense to look at their respective Fourier transforms/spectra.

For stochastic processes on the other hand, I am supposed to work with power spectral density in terms of qualitative analysis.

Why?

Because a stochastic process itself doesn't have a Fourier transform.

That's really all there is to it.

You can only transform signals (i.e. functions over a body isomorphic to $\mathbb R$, for example, functions of time). You can't transform a random variable whose individual realizations are such functions!

As @Marcus Müller mentions a stochastic process does not have a Fourier transform.

From a practical point of view, assume you have $N$ realizations of a signal, either deterministic or stochastic. Then, in order to estimate the power spectral density, it is common practice to average the power spectral densities of the individual realizations. For both deterministic and stochastic signals this will result in a more accurate estimate of the power spectral density as $N$ increases.

On the other hand, if you try the same with the Fourier transforms of the individual realizations this will only work for deterministic signals. For stochastic signals the phase of the Fourier transform is random and averaging generally result in a spectrum converging to zero as $N$ increases.


How to construct Polyphase (Frank) codes in MATLAB I have written a piece of code in MATLAB to generate Frank codes - Pastebin link for pulse compression.

I know that Frank codes are stepped frequency approximations to Linear Frequency modulation(LFM).

LFM Parameters are obtained here:

PW = 100e-6;% Pulse Width
PCR = 100; % Pulse Compression Ratio
CPW = PW/PCR;% Compressed pulse width
Fm = 2/CPW; %Frequency sweep

Stepped frequency approximation to LFM:steppedFreq = 0:Fm/(N-1):Fm;
The phase matrix is generated in the line:Matrix = mod(Matrix*incPhi,360);
Frank code is generated inside the for loop: FCP = [FCP,exp(1j*(2*pift + Matrix(i,p)*pi/180))];
The autocorrelation in dB is calculated here:AutoCorrelationdB = 20*log10(abs(AutoCorrelation)/max(abs(AutoCorrelation)));

My question:

The Peak side lobe level obtained is only around $3.5\textrm{ dB}$. Am I following a correct method for generating Frank codes? The peak side lobe level is given by the formula 20*log10(1/(N*sin(pi/M))) as per the book Principles of Modern Radar: Basic Principles by M.A. Richards - Page 825 which is equal to $29\textrm{ dB}$ for my case($N = 100,M = 10$).

% created by MUSTAFA SAMI AHMED
% Email mustafa_sami87@yahoo.com
clc;
clear all;
M =100;
L=sqrt(M);
bv = 0:L-1;
Matrix = bv'*bv ;
anglout = radtodeg(pi);
% alfa= (4*pi)/(L);

% Matrix_new = alfa*Matrix;
incPhi = (2*pi)/L; %incremental phase change
Matrix = mod(Matrix*incPhi,2*pi);

b=reshape(Matrix,1,[]); % matrix 1D size 1xM
for p = 1:M
  d(p) = exp(1j*b(p)) ;
end

AutoCorrelation = xcorr(d);

AutoCorrelationdB =10*log10(abs(AutoCorrelation)/max(abs(AutoCorrelation)));
 figure(1);
 stem(b);
%   x=wden(AutoCorrelationdB,'heursure','h','mln',10,'sym10');
  %x=wiener2(AutoCorrelationdB, [3 80 ]);
  figure(2);
plot((AutoCorrelationdB));


How to pad a windowed sinc filter in the frequency domain I'm creating windowed sinc filters to apply them to certain signals that I'm dealing with. To design the filters, I'm using the approach described in the book "The Scientist and Engineer's Guide to DSP". Here's a brief resuming:

\[h[i] =
\left\{
\begin{matrix}
Kw(i)\frac{\sin(2\pi f_c(i - \frac{M}{2}))}{i - M/2} & 0\le i \le (M+1)\\
\\
0 & \text{otherwise}
\end{matrix}
\right.
\]

\[
w(i): \text{A window function}\\ K: \text{DC gain}\\ f_c: \text{The cut frequency expressed as a fraction of the signal sampling rate}
\]

After calculating the filter kernel, I obtain the frequency response of the filter by taking the FFT of h[i] using an algorithm that approach the filter length to the nearest power of two(to get better performance). The same procedure it's also applied to the input signal in order to get its frequency response.

To my understanding, the next step would be multiply the two frequency responses to get the filtered signal in the frequency domain, but the problem is that both signals are not the same size and I haven't found an explanation on how to deal with this problem. So people, how should I pad the filter signal to make this point-wise multiplication? Should I fill the end of the filter with zeros to match with the signal size or should I do it symmetrically putting zeros to the left and right of the filter response?

Thanks in advance!

You don't pad it in the frequency domain, you pad it in the time domain (i.e. before you calculate the FFT).  You can put the zeros either symmetrically or after the filter.  Either way works, it just results in a shift in the output.

You will have to pad both $h[i]$ and the signal $x[i]$ with zeros so that both their lengths are atleast $M+N+1$ where $N$ is the length of $x[i]$. You can pad more so that the length of the padded vectors are powers of $2$.

This has to be done to ensure that circular convolution behaves like linear convolution.


Find Stopband Edge Frequency I tried to find the stopband edge frequency.I know that the signal amplitude falls to -3 dB below the passband edge frequency.I examine a lot of DSP books so i don't any information about how to calculate stopband edge frequency.

In general there is no standard way to calculate the stop band edge, other than specifying a desired stop band attenuation and checking at which frequency this attenuation is first attained (for a low pass filter).

For filters with an equiripple behavior in the stop band (e.g. FIR filters designed by the Remez exchange algorithm, or Chebyshev II and elliptic IIR filters), there is a clear stop band edge, which is the lowest frequency (for low pass filters) at which the attenuation equals the minimum stop band attenuation (attained at the local maxima of the filter's magnitude response in the stop band).

However, your filter looks like it has a monotonic magnitude response, so the definition of the stop band edge is arbitrary.

edit: This is for passband edge frequency.
One Way the cutoff frequency can be calculated as:
consider H(jw) normalized filter transfer function
$$
 H(jw)=\frac{V_0(jw)}{V_i(jw)}
$$
Here the 3dB frequency is when H(jw) is 1/sqrt(2) time the max value.Then it is finiding the values omega which satisfies the below equation
$$
H(jw)=\frac{1}{\sqrt2}
$$
This we can see when we take logarithm of 1/sqrt(2) we get approximately -3dB


About frequency resolution The frequency resolution of DFT is

$\Delta f = \frac{1}{T_{0}} = \frac{1}{NT} = \frac{f_{s}}{N}$

where $f_{s}$ is sampling frequency, $T_{0}$ is sampling time, $N$ is number samples, and $T$ is the sampling interval.

Where does this definition come from? Espescially the $1/T_{0}$, why this equality?

The longer the total sample window (T0), the smaller the frequency difference between two sinusoids with exactly M periods and exactly M+1 periods (for integer M <= N/2) in that window, which is required for orthogonality between adjacent DFT bins.

Also, the bin spacing is not the resolution. The frequency resolution can be greater or less than the DFT bin spacing, depending on the S/N ratio and your resolution needs (separation or estimation).

The resolution comes from the Rayleigh resolution - it is the distance from the peak of the main lobe of the response to the first null or zero in the response. In the continuous Fourier Transform this is a sinc() function.

Lets examine the DFT of a rectangular window (note - an arbitrary frequency can be created by multiplying by a complex exponential and using DFT transform properties i.e. frequency shifting).

For the time being lets ignore the sampling rate and just use indices in the time and frequency domain. So if we have a rectangular pulse of length $K$ samples and take a DFT of length $N$ then if $X(m)$ denotes the DFT of the rectangular pulse $x(n)$. Then (see here)

$$X(m) =e^{j\phi(m)} \frac{\sin(\pi mK/N)}{\sin(\pi m/N)}$$

The peak of the response occurs when $m=0$. The first null in the response occurs when
$$\pi m K/N=\pi $$ or $$ m=N/K.$$ Now to convert this index into a frequency - we know that the DFT frequency bin spacing is given by $f_s/N$ where $f_s=1/T$ is the sampling frequency and $T$ is time domain sample spacing. This gives:
$$ f_{res} = m\frac{f_s}{N} =\frac{N}{K} \frac{f_s}{N}=\frac{f_s}{K}=\frac{1}{KT}=\frac{1}{T_0},$$
where $T_0$ is the time duration of the signal (as per your original post).

Note - Using $\Delta f$ to denote resolution is problematic because it is also used to denote the DFT bin spacing - which is not the same thing as resolution.

Note - The phase term $\phi(m)$ is kind of complicated (see the link I provided), but it is irrelevant to this analysis because we are only interested in the magnitude of the response i.e. where is the peak and where is the first zero/null.


Principal Component Analysis as a step between STFT and MFCC For computing Mel Frequency Cepstral Coefficients you can use already calculated STFTs as a basis and perform the Mel frequency mapping on it.

My question: Does it make a significant difference if I calculate the STFTs, perform a PCA-transform on them and then calculate the MFCC compared to computing the PCA-transform at the very end on the MFCCs?

In both cases we have two dimensionality reductions (PCA and MFCC) and I'm not sure weather the order makes a difference.

One dimensionality reduction method (PCA) is linear, the other (MFCC) is  not. It  is thus unlikely that there is no difference at all since they probably do not commute. One serves more decorrelation, the other feature extraction.

The importance of the order thus depends a lot on your signal class, and what your purpose is.


Are reconstruction filter always needed? in mixed (digital and analogue application) where you have to convert a signal from continuous to digital time domain and then back again to continuous time an anti-aliasing LPF is needed and in general also a LP reconstruction filter.

Are there some applications (if so what are) in which the reconstruction filter is not needed or not strictly necessary?

If the system naturally has only a low enough frequency response (say a very heavy object relative to the forces involved), and the noise tolerance is high enough, then the effects of any aliasing, due to the lack of any additional low-pass and anti-aliasing filtering, may already be rolled-off below your desired or required noise floor.

No it not always needed.
For example in sound application that some weird sound is acceptable (not intended at DSP level).
Lets say I make up some sound using 8-bit microcontroller (without dedicated DSP module or DAC). I can just use R-2R Ladder to generate analog output from normal digital output pin. The output will not smooth at all but it is acceptable (It's art).
Depend on application requirement sometimes you may not need it.


Fourier transform of discrete time unit step function To obtain fourier transform of u[n],
u[n] - u[n-1] = delta[n] , taking fourier transform of both sides of the equation results in :
U(w) - exp(-jw) U(w) = 1 , hence :
U(w) = 1/(1-exp(-jw)) which is wrong and the right answer has an extra term.
Which step is wrong in this possible solution?
I know the right proof of fourier transform of u[n], my question is regarding the wrong part of this solution.

The DFT of a unit step response is
$$U(\omega) = \frac{1}{1 - e^{-j \omega}} + \pi \delta(\omega)$$
Applying the shift property as you did will give:
$$\mathcal{F}(u[n] - u[n-1]) = U(\omega) - U(\omega)e^{-j \omega} = \frac{1}{1 - e^{-j \omega}} + \pi \delta(\omega) - [\frac{1}{1 - e^{-j \omega}} + \pi \delta(\omega)]e^{-j \omega}$$
that is
$$\mathcal{F}(u[n] - u[n-1]) =  \frac{1- e^{-j \omega}}{1 - e^{-j \omega}} +\pi\delta(\omega)( 1 - e^{-j\omega}) = 1+\pi\delta(\omega)( 1 - e^{-j\omega})$$
The second term is always zero because for $\omega = 0$, $1 - e^{-j\omega} = 0 $ and it's zero on any other point, so you get,
$$\mathcal{F}(u[n] - u[n-1]) = 1 = \mathcal{F}(\delta[n])$$


Image Projection with SIFT Recently, I am thinking about a problem regarding image projection. I have a larger picture as bellow.

I want to find the transformation matrix, x'=Hx. When the above image pass the transformation matrix H, I will get another smaller picture as bellow.

I have tried to use SIFT to find four major points and use these points to implement DLT algorithm. However, after executing SIFT, I find that there are many matching error as bellow. The distRatio is 0.5

I have thought for almost two weeks. Is there any different ideas which can be use to solve this problem. Thank you very much.

The logo you are looking at is symmetric. For example, you could rotate it by 120 degrees and get the same pattern. If you take a look, many of the "errors" are not really errors. The same pattern appears in your image many times. It is called self-similarity.

You could do several things:

Change descriptors to avoid confusion between features that are similar in an image.
Instead of finding one match for keypoint, find 3 matches and test multiple combinations by fitting a transformation and checking which fits better.
Limit the transformation by rotation angle, mirror, etc. For example, angle must not be larger than 60 degrees


Upscaling detection approach I'm trying to implement a solution that would, to some extent, evaluate a picture telling what are the chances that it was upscaled from the smaller image.

My approach is as follows:

Take a picture of unknown origin (scaled or not)
Perform sharpening using some technique which sharpens only evident edges (1)
Extract edges using Difference of Gaussians
Measure entropy of resulting image

Then downscale the image and upscale it back to the same size and compare entropies.
OR
Then upscale the image and downscale it back to the same size and compare entropies.

Does it sound reasonable? What technique should I use as (1)? Also, maybe I should do it in small steps so I don't introduce too many artificial points? I mean if the image was upscaled (preferably instantly, to target size) it won't matter a lot - but it the image weren't upscaled this incremental improvement will be better for the picture than taking it and upscaling instantly?

Here is the poor man's approach:

var histogram = []
foreach pixel
    histogram.append(abs(pixel.value - pixel.south.value))
    histogram.append(abs(pixel.value - pixel.east.value))
sort(histogram)
sharpness = histogram[histogram.count * 0.95]
upscaled = sharpness < X // You choose X experimentally

A simple method for comparing images is template matching

Its usually used for finding a smaller image in a larger image.

You could also use it as a way of comparing how similar two images of the same size are.

The opencv article I linked to lists a number of different correlation functions such as sum of squared differences and correlation coefficient.

Whether or not this will be accurate enough for what you are doing will only become clear after you have tried it.


Digital Image Processing I am trying to extract only urban aeas from satellite imagery Using Matlab. I have applied some edge detectors like canny and also used morphological functions and I am able to extract all the boundaries ,but not able to seperate them from various urban aeas like roads, trees and waterbodies.

Please help i can share the code if any one wants

I think this code would be helpful, as displaying intensity part gives you the spatial characteristics of an image.

%Write a MATLAB function which converts an RGB to HSI. Display the Hue
%image, Saturation image and the Intensity image. Histogram equalize the H,
%S and I images and the HSI image. Then convert the HSI image back to RGB.

function rgbtohsi(x)
F=imread('C:\Documents and Settings\admin\Desktop\TM\tanya\img3.tif ');
F=im2double(F);
r=F(:,:,1);
g=F(:,:,2);
b=F(:,:,3);
th=acos((0.5*((r-g)+(r-b)))./((sqrt((r-g).^2+(r-b).*(g-b)))+eps));
H=th;
H(b>g)=2*pi-H(b>g);
H=H/(2*pi);
S=1-3.*(min(min(r,g),b))./(r+g+b+eps);
I=(r+g+b)/3;
hsi=cat(3,H,S,I);
HE=H*2*pi;
HE=histeq(HE);
HE=HE/(2*pi);
SE=histeq(S);
IE=histeq(I);
choice=input('1:RGB to HSI\n2:Display Hue, Saturation and Intensity Images\n3:HSI to
    RGB\n4:Hue-Equalization\n5:Saturation-Equalization\n6:Intensity-
    Equalization\n7:HSI-   Equalization\n Enter your choice :');
switch choice
case 1
    figure,imshow(F),title('RGB Image');
    figure, imshow(hsi),title('HSI Image');
case 2
    figure,imshow(F),title('RGB Image');
    figure, imshow(H),title('Hue Image');
    figure, imshow(S),title('Saturation Image');
    figure, imshow(I),title('Intensity Image');
case 3
    C=hsitorgb(hsi);
    figure,imshow(hsi),title('HSI Image');
    figure, imshow(C),title('RGB Image');
case 4
    RV=cat(3,HE,S,I);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-Hue Equalized');
case 5
    RV=cat(3,H,SE,I);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-Saturation Equalized');
case 6
    RV=cat(3,H,S,IE);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-Intensity Equalized');
case 7
    RV=cat(3,HE,SE,IE);
    C=hsitorgb(RV);
    figure,imshow(hsi),title('HSI Image');
    figure,imshow(F),title('RGB Image');
    figure, imshow(C),title('RGB Image-HSI Equalized');
    F1=C;
    otherwise
    display('Wrong choice');
end
end
    function C=hsitorgb(hsi)
    HV=hsi(:,:,1)*2*pi;
    SV=hsi(:,:,2);
    IV=hsi(:,:,3);
    R=zeros(size(HV));
    G=zeros(size(HV));
    B=zeros(size(HV));
%RG Sector
    id=find((0<=HV)& (HV<2*pi/3));
    B(id)=IV(id).*(1-SV(id));
    R(id)=IV(id).*(1+SV(id).*cos(HV(id))./cos(pi/3-HV(id)));
    G(id)=3*IV(id)-(R(id)+B(id));

%BG Sector
    id=find((2*pi/3<=HV)& (HV<4*pi/3));
    R(id)=IV(id).*(1-SV(id));
    G(id)=IV(id).*(1+SV(id).*cos(HV(id)-2*pi/3)./cos(pi-HV(id)));
    B(id)=3*IV(id)-(R(id)+G(id));
%BR Sector
    id=find((4*pi/3<=HV)& (HV<2*pi));
    G(id)=IV(id).*(1-SV(id));
    B(id)=IV(id).*(1+SV(id).*cos(HV(id)-4*pi/3)./cos(5*pi/3-HV(id)));
    R(id)=3*IV(id)-(G(id)+B(id));
    C=cat(3,R,G,B);
    C=max(min(C,1),0);
end


parametric eq with FIR filters Are there any good published approach for designing parametric EQs as FIR filters.  The standard techniques discussed in DSP textbooks are frequency-domain-sampling or using Parks-McClellan, which this involves passband, stopband, transition band, etc.  This is nice for rejecting a certain portion of the frequency spectrum, or mimicking and existing analog/digital filter.  However what if I want to design a filter directly as an FIR filter, for example a peaking/notching/shelving/low-pass/high-pass filter which takes frequency, gain, and bandwidth parameters?  Is there a standard technique for designing parametric EQs directly as FIR filters?

These filter types are well covered here http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

However these are all single biquad IIR implementations. To turn this into an FIR, you could simply calculate the impulse response of these filters and window to a finite size with the desired amount of accuracy. Another option is to design the filter as IIR, sample the frequency response and then do any of the many FIR fitting techniques.

In general you will find that FIR filters are not well suited for this and the number of taps required will heavily depend on the relationship between the corner frequency and the sample rate and also on the Q of the filter.

Matlab's documentation for filter design should have the references you are looking for.

For example fir2 might do it.


Find the $\mathcal Z$-transform of this function? I need to find the $\mathcal Z$-transform of $x(n)=a^{-n} u(n)$. Assume, $a$ is a positive constant , but the power of $a$ is negative.

Looking at the transform table, I found that $\mathcal Z$-transform for $x(n)=a^{-n} u(n)$ is not there, so I tried to found it by on my own, but I am not sure if my answer is right, please can anyone check it and correct me, in case?

Doing some calculations is appreciated. However, it is preferred that you write them in the question because a picture can vanish more easily. Do not hesitate to edit the math.

To help you with an hint in checking your answer, remember that:

$$ a^{-n} = \left(\frac{1}{a}\right)^n$$

a standard trick to play with known tables.

No you have done a small mistake how a^-n became a^n? So it is like this the ROC is right sided


What are poles and zeros? The concept of poles and zeros in filters was introduced to me during our lab (our lab isn't sync with our lecture) through the pole-zero plot generation of filters in MATLAB. I didn't get its importance until it was mentioned in our lecture that poles make the equation infinity (or its denominator zero) while the zeros make the roots in the numerator zero. I still can't quite grasp what poles and zeros are. Can someone give an an intuitive and simple explanation on what poles and zeroes are? Much appreciated.

Take the equation b/(x-c) with b non-zero.  The ratio goes to infinity as x approaches c.  So c is the location of a pole (something tall and pointy in a graph).

Take the equation (x-b)/c with c non-zero.  The ratio goes to zero as x gets closer to b.  So b is the location of something commonly called "a zero".

You can not only do this with scalar x, but with complex x, thus the domain of the poles and zeros will be on a (complex) plane instead of on a line.

If the ratio represents something about the response of a filter, it might say the filters output is at or near zero when the input is at or near a response "zero".  And bad things might start to happen when x gets close to a pole (power supply starts smoking when asked to supply infinite amps, math operations produce NaNs or fixed-point overflows, etc.)

First, you should see a z-plane as a set of complex exponential signals. If $z = z_0$, then corresponding discrete signal $Re(z^n) = Re(z_0^0),\ Re(z_0^1),\ Re(z_0^2),\ ...$. If $|z_0|<1$, it's a decaying signal.

Second, poles and zeros are used to describe an IIR system, i.e. a system with a feedback.

Zeros are easy: if system has a zero at $z_0$, that means a signal defined by $z_0$ on a z-plane will pass through a feedback loop and sum with itself strictly out-of-phase resulting in zero output. Poles are a bit trickier: if system has a pole at $z_0$, that means a system will generate this signal than it's disturbed and moving freely. Obviously if system has a pole for increasing signal ($|z_0|>1$), it's unstable. When designing a filter, zeros should be placed to suppress undesired frequencies. Poles are placed such that

they cancel out the the influence of zeros for the frequencies to be passed
filter remains stable


Checking if u[n] and v[n] are orthogonal Let $x[n]$ and $y[n]$ be two orthogonal baseband signals, with bandwidth B. Check if the following signals are orthogonal:

$u[n] = x[n]*\cos(2*\pi*f_c*n)$

$v[n] = y[n]*\cos(2*\pi*f_c*n)$

I'm not sure what the problem asks to be honest.
1) What exactly does he mean by bandwidth?? I only know the definition of a bandwidth of a sinc function, but how can I think of it generally?

2) I assumed the case for a rectangular function and reduced the problem to:

$\sum_{n=-\frac{1}{B}}^{n=\frac{1}{B}}\frac{1}{2}\cdot x[n]\cdot y[n]\cdot \cos(4\cdot \pi\cdot f_c\cdot n)$

But I don't know where to go from there. It's also weird for me to evaluate since $\frac{1}{B}$ is likely not an integer.

EDIT: also, $f_c$ >> B

It's essentially asking you to show that two signals that are orthogonal in the baseband (their original form) remain orthogonal after amplitude modulation process. In other words, amplitude modulation preserves orthogonality.

Think about the question this way. You have two signals that are orthogonal, meaning that

$$ \sum_{m=-\infty}^{\infty}x[m]y[m] = 0 $$

This is a bit ugly, since we have to deal with infinities. However, your problem also says that they both have bandwidth $B$, which means that in frequency domain they will have finite width.

Let's define $\mathscr{B}$ as

$$ \mathscr{B} = \frac{B}{F_s}$$

which means that in your discrete Fourier transform signals $X[k]$ and $Y[k]$, the only non-zero samples exist for $-\frac{\mathscr{B}}{2} \le k \le \frac{\mathscr{B}}{2}$. Since DFT is an orthogonal transform, signals orthogonal in time domain will remain orthogonal in frequency domain, i.e.

$$ \sum_{l=-\infty}^{\infty}X[l]Y[l] =\sum_{l=-\mathscr{B}/2}^{\mathscr{B}/2}X[l]Y[l] = \sum_{m=-\infty}^{\infty}x[m]y[m] = 0 $$

We know that multiplying a signal by a cosine shifts is away from DC in frequency domain, and the fact that the modulating frequency $f_c$ is larger than $B$ means that this won't cause any problems with overlapping modulated bands (think about what would happen if this weren't the case).

I think this gives you enough hits to take it from here.

Let x[n] and y[n] be two orthogonal baseband signals, with bandwidth B. Check if the following signals are orthogonal:

$u[n]=x[n]∗cos(2∗π∗fc∗n)$
$v[n]=y[n]∗cos(2∗π∗fc∗n)$

Answer:
For $u[n]$ to be orthogonal to $v[n]$
Inner Product = $<u[n],v[m]>$

By using Parsevals theorem
Compute
Integration (from -Inf to Inf) {U(f) * conj(V(f))} df

Here $U(f)$ is frequency domain of $u[n]$
= $(X[f-fc]+X[f+fc])/2$

Here $V(f)$ is frequency domain of $v[n]$
= $(Y[f-fc]+Y[f+fc])/2$

Inner Product will be
= Integration (X[f-fc].Y[f-fc] + X[f+fc].Y[f+fc])/4
= 0


Difference between Digital signal and Discrete signal I have just started studying Digital Signal Processing. Can someone explain what is the difference between a Discrete Signal and a Digital Signal in simple words?

Thanks in advance !

Suppose you have a continuous time analog signal. It is continuous in both time and amplitude. Now when you sample it ,you get discrete samples every Ts seconds. Now you have discrete samples(discrete in time) each of which can take continuous value( in amplitude). This is normally referred to as discrete signal(discrete in time but continuous in amplitude).

Now further when you take this discrete signal and quantise it i.e. say assign each sample which takes continuous amplitude to one of the N discrete quantisation levels of a quantiser, Then your total signal now is a digital signal. So a digital signal is discrete in time and discrete in amplitude.

Digital signal is a continuous signal.Discrete signal is a non continuous signal.
i.e the digital signal is present at all time(continuous signal). the discrete signal is present some time values only(non continuous signal)


How do you relate imaginary numbers with phase shift? How to imagine this? So I found this answer on one of the questions. I wanted to ask how the complex part can be related with the phase shift. I wanted to ask this on the same question but it didn't allow me. Could anyone please help me understand this concept. Even if I imagine this as a spiral with the value of the imaginary part changing, the starting time of the wave remains the same doesn't it? So how is it related to phase shift?

"Negative frequency doesn't make much sense for sinusoids, but the Fourier transform doesn't break up a signal into sinusoids, it breaks it up into complex exponentials:

$F(\omega)=\int_{-\infty}^{\infty} f(t)e^{−jωt}dt$

These are actually spirals, spinning around in the complex plane:

complex exponential showing time and real and imaginary axes

Spirals can be either left-handed or right-handed (rotating clockwise or counterclockwise), which is where the concept of negative frequency comes from. You can also think of it as the phase angle going forward or backward in time.

In the case of real signals, there are always two equal-amplitude complex exponentials, rotating in opposite directions, so that their real parts combine and imaginary parts cancel out, leaving only a real sinusoid as the result. This is why the spectrum of a sine wave always has 2 spikes, one positive frequency and one negative. Depending on the phase of the two spirals, they could cancel out leaving a purely real sine wave or a real cosine wave or a purely imaginary sine wave, etc."

Paint a red dot on the blue complex spiral at t=0, where the value is 1.0.  Now roll the spiral upside-down until the red dot is at -1.0, still at t=0.  This roll would be equivalent to a phase shift of Pi.  Note that the real part no longer corresponds to a cosine wave that starts at t=0.

But the real part of the rolled spiral does correspond to a cosine wave that starts (with value 1.0) at +Pi or -Pi, which is a different starting time, contradicting one of your assumptions above.

Thus phase is related to starting time.

I think the 3D plot is making this concept more difficult than it really is.

If your signal is f(t) = cos(2πft), then the phase of f(t) is simply 2πft. The real part of the vector is cos(2πft) and the imaginary part is sin(2πft).

I am not sure why someone would use a 3D plot to describe this. These concepts remain on the complex plane, and the angle of the vector (its phase) is simply 2πft.

The 3D plot you show is how we explain the phase of a propagating plane wave. Then the distance between each loop of the spiral is called the wavelength and is equal to c*T where c is the speed of light and T is the period of the signal. These concepts are important to RF engineers for example, who need to understand the phase of the voltage and current at each point in space as the signal propagates.

I think you will have less difficulty if you stay on the complex plane.


Is there a difference between Image Alignment and Stereo Rectification for stereo correspondence or matching? As far as I know,
Image alignment just finds the best fit given two images of the same scene.

Whereas Image rectification, with regards to stereo vision, warps the images so that the epipolar lines are all on the same level. Essentially all common pixels between the two images are warped to be on the same y axis.

My question is, what is the difference between Image Alignment and Image Rectification, given a pair of images, when doing stereo correspondence?

If I were to find the disparity between two images, how is image alignment different to image rectification?

There two ways to look at this problem.

In simple terms, image rectification warps both images onto a common coordinate frame by typically estimating the transformation using the epipolar geometry.
Image alignment finds the transformation from one image to the other. It doesn't guarantee any constrains on the epipolar geometry and only one single image is warped.

Image rectification might benefit from image alignment (of local patches) to obtain the correspondences. Once these are found, the epipolar lines can be forced to be parallel.

Well image alignment is quite general, some techniques for get the disparity use image alignment in small regions, and they can estimate the disparity by using those methods. And depending on the technique of image alignment is what you are going to obtain, if you are using the borders of the image then the process is going to find the best match for the edges no matters if the images looks very similar or not. And you have to consider that in stereo pairs you are going to have occlusions in images that can be translated to missmatching issues.