Hilbert transform properties Here Its says Hilbert transform is a non-causal,linear,and time-invariant system

How can I prove it mathematically?

wikipedia says the input output relation like this $$\boxed{y(t)=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(\tau)}{t-\tau}d\tau}$$

so from this relation it showing time varying nature because

for $X(t-t_o)$,$y(t)$ is $${y(t)=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(\tau-t_o)}{t-\tau}d\tau}$$

and $${y(t-t_o)=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(\tau)}{(t-t_o)-\tau}d\tau}$$

So both are not same so its Time variant

HINT:

Prove that
$$y(t)=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(\tau)}{t-\tau}d\tau=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(t-\tau)}{\tau}d\tau\tag{1}$$
Use the right-most expression in $(1)$ to show that the system is time-invariant.

Equivalently, you can directly show that the original input-output relation is a convolution integral, from which it follows that the system must be LTI. Then find the expression for the impulse response $h(t)$ and show that $h(t)\neq 0$ for $t<0$, which means that the system is not causal.

EDIT:

The same variable substitution that can be used to show $(1)$ can also be used to directly show the equivalence of the response to $x(t-t_0)$ and the delayed response to $x(t)$:

$$y_2(t)=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(\tau-t_o)}{t-\tau}d\tau$$

$$\begin{align}y(t-t_o)&=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(\tau)}{(t-t_o)-\tau}d\tau{\Huge|}_{\tau=u-t_0}\\&=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(u-t_0)}{(t-t_o)-(u-t_0)}du\\&=\frac{1}{\pi}\int_{-\infty}^{+\infty}\frac{x(u-t_0)}{t-u}du\\&=y_2(t)\end{align}$$


Median Filter for salt and pepper noise removal Why median filter is considered as good for removal of salt and pepper noise? What are the other filters used for the same?

Median filter is considered good because unlike averaging filter which ruins the edges of an image by blurring it to remove the noise, median filter removes only the noise without disturbing the edges.

Well, median filter is the best and only filter to remove salt and pepper noise.

Hope this helps:)

Thank you!!!


Designing ideal lowpass filter in the form of FIR How to design a ideal lowpass filter by cutoff frequency of π/4 in the form of FIR filter by 64 samples using matlab?

You  don't.

An ideal low pass FIR filter has infinite length so the requirements of "ideal low pass" and "64 taps" are mutually exclusive.

You can approximate an ideal filter, but the best way to do this depends on the specific requirement and trade-offs of your application. This being said, Richard's answer is a really good starting point :-)


Is the common FIR Filter delay line (see pic) also valid on purely logical FPGA implementations? Lets say we want implement simply the mean with constant coefficients. So all b's are set to e.g. b=0.5. In this case I do not understand the delay line in the picture below.

I mean I understand it for von Neumann, Harvard architectures.., but how about if I just calculate everything on parallel FPGA? Shouldn't be possible to implement the whole structure without any delays? (..because it is all parallel and there is no feedback) or is it in practice just not possible?

http://en.wikipedia.org/wiki/Finite_impulse_response#/media/File:FIR_Filter.svg

(What I would understand is that there is maybe a delay at the end when the summation takes place)

I'm not sure if I understood our question right, and therefore sorry if I point something wrong, but the delay in the delay line comes from signal sampling and not from computational processing type. Therefore whether you realize that signal flow graph in a von Neumann type serial architecture or in aparallel FPGA based architecture doesn't matter. You will always have the delay per input sample as dictated by your systems sampling rate.

On the other hand of course a parallel implementation would take much less clock cycles, based on the amount of parallelisation, to compute per sample output campared to a serial architecture of a general purpose CPU, however as I said the delay in the delay line is nothing related with computational delay of the hardware system.

Finally the signal flow graph yields a pictorial representation of processing structure, you don't really implment those delays. They are represented by registers or memory locations, as possibly you know, which are updated either per sampling clock or via buffer copies which provides that delay effect represented in the flow graph.


Deriving the integration property of the Fourier Transform I want to derive the property of the Fourier Transform that states that if $X(j\omega) = \mathcal{F} (x(t))$ then $$\mathcal{F} \left( \int_{-\infty}^{t} x(\tau) \mathrm{d} \tau \right) = \frac{1}{j\omega} X(j\omega) + \pi X(0) \delta(\omega).$$
I started as follows, since $$x(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega) e^{j\omega t} \mathrm{d}\omega,$$ integrating both sides gives us $$\int_{-\infty}^{t} x(\tau) \mathrm{d} \tau = \int_{-\infty}^{t} \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega) e^{j\omega\tau} \mathrm{d}\omega\mathrm{d} \tau$$ $$\int_{-\infty}^{t} x(\tau) \mathrm{d} \tau = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega) \left(\int_{-\infty}^{t} e^{i\omega \tau} \mathrm{d}\tau \right) \mathrm{d}\omega.$$

How do I evaluate the inner integral inside the parantheses? My gut feeling says that it's related to the Fourier Transform of the constant function $y(t) = 1$ but I can't relate it to that since the upper limit is $t$ rather than $\infty$. How do I reach the result above?

Thank you in advance.

First, recognise the integral of $x(t)$
$$\int_{-\infty}^{t} x(\tau) d\tau \tag{1}$$
as a convolution with $u(t)$
$$ \int_{-\infty}^{t} x(\tau) d\tau =x(t) \star u(t) \tag{2}$$
where $u(t)$ is the unit-step function.
Then from the convolution property of the CTFT, one gets the following :
$$ \mathcal{F}\{ \int_{-\infty}^{t} x(\tau) d\tau  \} = \mathcal{F}\{ x(t)\star u(t) \} = X(\omega)U(\omega) \tag{3}$$ where $U(\omega)$ is the Fourier transform of the unit-step function which is $$U(\omega) = \frac{1}{j\omega} + \pi \delta(\omega) \tag{4}$$
Then the Fourier transform of the integral of $x(t)$ is
$$ \mathcal{F}\{ \int_{-\infty}^{t} x(\tau) d\tau \} = \frac{ X(\omega) }{j\omega} + \pi X(0) \delta(\omega) \tag{5}$$


Calculate the power of a discrete-time signal We got this question in our test :-
$$S(t) = 5\cos(10πt - π/2) + 6\sin(15πt)\cos(15πt)$$
My Solution :-
$$S(t) = 5\cos(10πt - π/2) + 3\sin(30πt)\qquad    \textrm{[Using $\sin(2A) = 2\sin(A)\cos(A)$]}$$
=> Power$ = 5x5/2  +  3x3/2
         = 25/2 + 9/2
         = 17$
But in the solution, the correct answer was given to be 21.5 watt.
What mistake did I do?

Apart from possible dimensional problems as discussed in the comments, your result is correct and the given solution is wrong. In general with such exercises, don't forget to check if the cross-term really cancels out (in this case it does).


Conditions for expressing a system as difference equation I have read when studying this subject that a system is LTI and causal if and only if it can be expressed as a difference equation (if it is in continuos time, as a differential one). I don't know if this is true, mayne I'm remembering it wrong.

The problem is that today I thought of the case of an IIR filter. It is not causal and it can be written as a difference equation.

So I have two questions about this:

1) Why does an LTI system have to accomplish the initial rest condition (causality?) in order to be expressed as a difference equation?

2) What about the IIR filters? They are not causal and they can be expressed as a difference equation. Why is that? What am I thinking wrong?

A difference equation does not imply causality.

For example, consider this system, with input $x$ and output $y$:

$y[n] + 0.2 y[n-1] = x[n] + 0.5 x[n-1]$

You can solve it (recursively) forwards in time as $y[n] = x[n] + 0.5 x[n-1] - 0.2 y[n-1]$.

Or you can solve it backwards in time as $y[n-1] = 5 x[n] + 2.5 x[n-1] - 5 y[n]$.

These are 2 valid and different results for the same difference equation. These 2 results correspond to the 2 possible regions of convergence of the system response (Z transform). Inner ROC corresponds to the causal solution; outer ROC to the non-causal solution.

Systems of greater order will have more ROCs and more different solutions.

If you impose causality, then you are implicitly selecting the inner ROC.

For the second question, IIR filters can be causal. The example above is an IIR system. The causal impulse response will be $h_{causal}[n] = (-0.2)^n u[n] + 0.5 (-0.2)^{n-1} u[n-1]$. This is causal and IIR.


Conceptual Question on equalization technique in rayleigh fading channel based on a paper Question   Link http://www.raymaps.com/index.php/theoretical-ber-of-m-qam-in-rayleigh-fading/

gives theoretical expressions for BER in Rayleigh fading channel. Can I  use the expression for 64-QAM in this paper? The Authors have given a plot of the BER performance but there is no expression for it. So, I am unsure how to get the plot for Rayleigh case.

The model they have used is general and you can apply it to QAM as well.
In wireless communication systems, channel and noise are two different impairments affecting the overall transmission. Here as everywhere else, the noise is AWGN (additive white Gaussian; if you want to know what these terms mean, see this article). And the channel is a linear filter with taps whose magnitudes are modeled as Rayleigh random variables.
Yes, you can use the expression in that paper. The authors didn't give the expression because they might not have known it and simply generated it from Monte Carlo simulations.


Is it valid to calculate the magnitude, power, and phase of a real time-domain signal without converting it to the frequency-domain? I would assume you perform the calculations the same way, but since there would be no imaginary component because the signal is not complex it would be simpler:

magnitude = |x(n)|

power = |x(n)|^2

phase = 0 (for all values of n) - since there's no imaginary portion the phase equation would always result in a 0 value.

I'm trying to fully grasp and internalize some of the basic concepts of signal processing, and while I haven't ever really seen anyone try to use the equation for phase on a time-domain signal I don't see why it wouldn't be valid even if the results don't have much meaning.

No, this doesn't make much sense. What you can do in the time domain is compute the analytic signal and derive the signal's instantaneous amplitude (envelope) and its instantaneous phase from it.

Take as an example

$$x[n]=A\sin(\omega_0n+\theta)\tag{1}$$

The corresponding analytic signal is

$$x_a[n]=-jAe^{j(\omega_0n+\theta)}\tag{1}$$

Its instantaneous amplitude (envelope) is $\big|x_a[n]\big|=A$, and its instantaneous phase is $\arg\{x_a[n]\}=\omega_0n+\theta-\pi/2$.


Using Least Squares for curve fitting in Matlab I have a problem right now to solve. I have a Gaussian peak, which I need to divide into three components (three Gaussian distribution). For those three component Gaussian distribution, I already know the mean and standard deviation, if I want to do this by using least squares, where should I start with in Matlab?

Thanks

1) Form a matrix, A of order $m$ x $3$, such that:

$$A=\left[f(x) \ \ g(x) \ \ h(x)\right]$$

where $f(x), g(x)$ and $h(x)$ are three column vectors such that each column is a sample of component guassian at $x$.

2) Form a column vector, $b$ of your original Gaussian pulse of order m.

3) Solve the system $x=((A'*A)^{-1}*A')*b$, where $x$ is your desired set of coefficients of order $3$ x $1$, that would minimize the error in the least square sense.


Analytically determine a PSD from a transient function This question is related to a series of questions I have asked about the units of PSD and ESDs. I include it as a separate question as it may have worth in isolation.

As I understand it to compute the power spectrum density from some transient signal one takes the approach of computing the truncated Fourier transform
$$\mathcal{F}_T(\omega) = \int_{-T/2}^{+T/2} e^{-i \omega t} f(t) \ \ dt$$
Taking the absolute value square and then taking the limit as $T\rightarrow \infty$
$$\text{PSD} = S_{x}(\omega) = \lim_{T\rightarrow\infty}\frac{1}{T}|\mathcal{F}_T(\omega)|^{2}$$
When ignoring the normal dimensioning of units as described here by a very nice answer by Matt L., and also descibed here and in the Wikipedia article here, we get the familiar units of a PSD of $\rm{[V^{2}/Hz]}$.

My question

So that is the background and I believe I follow it (please say if there is a mistake or something overtly wrong!)

But If I want to analytically derive a function for a PSD or $S_{x}(\omega)$, I always end up with zero as the $\lim_{T\rightarrow\infty}$ kills the expression because of the factor of $1/T$. Is this the correct approach to analytically deriving a function or lineshape for a PSD with correct units? The end goal is to use it to fit a function to spectra acquired with an FFT spectrum analyser.

Energy signals, i.e., signals with finite energy

$$\int_{-\infty}^{\infty}|f(t)|^2dt<\infty$$

have zero power and, consequently, a power spectrum that is equal to zero. They do have an energy density spectrum, which is the squared magnitude of their Fourier transform.

You will only get a non-zero power spectrum, according to the definition in your question, for power signals, which have a finite non-zero power defined by

$$\overline{f^2(t)}=\lim_{T\to\infty}\frac{1}{T}\int_{-T/2}^{T/2}|f(t)|^2dt$$

I.e., power signals satisfy

$$0<\overline{f^2(t)}<\infty$$


How we can use the Hadamard transform in feature extraction from an image? Which coefficients of the Hadamard transform contain the highest energy compaction? Is there any regular pattern in them as for the  DCT?

The Hadamard transform is a linear and orthogonal transformation. It belongs to a wider class, under the names of Paley, Walsh, Rademacher and Hadamard. The ordering of the rows or columns of an Hadamard matrix is not well suited to a frequencial interpretation, notably because of the standard recursive definition (for Hadamard matrices of order $2^m$):

$$H_m = \frac{1}{\sqrt2} \begin{pmatrix} H_{m-1} & H_{m-1} \\ H_{m-1} & -H_{m-1} \end{pmatrix}$$

It is illustrated in the center matrix of the picture below.
However, the Walsh ordering (right most matrix) order the rows by the number of sign changes in the sequence (transition between the green and the red): the first row is of constant sign, the second one has only one sign  change, etc. It can be interpreted as related to the sign of a basic vector of the DCT (roughly), and then has "regular patterns". One then talks about the notion of sequency instead of frequency.

Basically, Hadamard and Walsh have the same compaction properties, in the sense of the $\ell_0$ count: for each coefficient with some energy in Hadamard, there is a coefficient with the same energy in Walsh.

However, for real data, the high energy coefficients tend to be scattered for Hadamard, and better localized toward some sequency band  for Walsh. For low-sequency signals, Walsh can be equipped with a zig-zag like ordering.

The quality of compaction depends on the data. For checkerboard-like patterns, Walsh-Hadamard can be better than the DCT. For a Gaussian white noise, both would perform equally (as for any orthogonal transform). For correlated signals, esp. whose autocorrelation can be modelled by an simple autoregressive model, the DCT tends to perform better. The relative performance of the DCT wrt Haar, Fourier, Walsh-Hadamard and Karhunen-Loève transforms is assessed for instance since Discrete Cosine Transfom, Ahmed et al., 1974. The DCT enjoys a close to optimal (Karhunen-Loève) performance, while retaining fast algorithm. You could find an adapted ordering for Hadamard, to keep the high-energy coefficients, but using Walsh is more convenient, as coefficients will tend to group in blocks, more than with Hadamard.

Finally, when even faster is needed, especially when the data is mildly correlated, Walsh or close-to-Walsh transforms (requiring no multiply) are quite effective (for instance in the residue after block prediction in video compression).


Periodicity of the output of causal time-invariant systems when the input is periodic For time-invariant systems, if the input is periodic, the we can argue that the output is periodic:

If we shift the input signal one period, we come up with the same signal as before, and since the output must be shifted the same amount as a result of the system being TI, it should be periodic too.

For causal systems we know that the output signal at time $t$ depends on the input for times up to $t$. My question is, can we conclude that for causal TI systems if the input is periodic before $t_0$, the output will be periodic up to $t_0$ too?

Your line of reasoning is correct. For a (not necessarily linear) time-invariant system the following always holds: if $y(t)$ is the response to an input signal $x(t)$, then $y(t+T)$ must be the response to the shifted input signal $x(t+T)$. If the input is periodic with period $T$, this implies that $x(t)=x(t+T)$, and, consequently, $y(t)=y(t+T)$, i.e. the output signal is periodic with the same period as the input signal.

Furthermore, if the system is causal and the input up to time $t_0$ is periodic, the output up to time $t_0$ must be periodic too, because a causal system cannot know the (possibly non-periodic) input beyond time $t_0$.

This is (part of the) theory. But then comes @MBaz with his "counter-example" of an ideal integrator. The problem here is indeed, as already noted in a comment, the integrator's instability. Showing the time-invariance of a linear system involves manipulating the convolution integral. However, if for the given input signal that integral is divergent, the proof becomes invalid. And this is the case with systems having poles on the $j\omega$-axis (or, for discrete-time systems, on the unit circle). If an LTI system has a pole on the $j\omega$-axis and if you excite it with the corresponding pole frequency, you get an output signal with increasing amplitude, which is obviously non-periodic. You can also see it in the following way: the response of an LTI system to a complex exponential $e^{j\omega_0 t}$ is

$$y(t)=H(\omega_0)e^{j\omega_0 t}\tag{1}$$

Eq. (1) of course requires that $H(\omega_0)$ can be evaluated. For an ideal integrator this is not the case for $\omega_0=0$. Note that if the input signal to the integrator had no DC component, then the response to a periodic input signal would indeed be periodic.


How can I generate time-variant Rayleigh channel in MATLAB One of known models for channel simulation is called Rayleigh Fading channel, it's usually used to simulate wireless channel in time variant/invariant.

As I know, Rayleigh fading channel can be considered as sum of two Gaussian random variables, we can use the below command to generate that channel:

h = 1/sqrt(2)*[randn(1,M+1) + 1i*randn(1,M+1)];

Where M is channel length and h will be the channel.

But what's about if I want to generate Rayleigh channel in time-variant environment?   is that possible?  Please correct me if the above command generates Rayleigh channel or there is another command.

thank you

Yes, you can have time variant channel in MATLAB either by using comm.RayleighChannel()  , you can read about it in help of matlab.

or the other way you use the same command you have mentioned but with each symbol the channel should be changed. so it will be like you will generate a channel which is Quasi-static. You above command will become:

h = 1/sqrt(2)*[randn(n_symbol,M+1) + 1i*randn(n_symbol,M+1)];

where n_symbol represents the number of symbols you have. then to do the convolution, must be use the loop of for.

for jj = 1:nSym

  new_signal(jj,:) = conv(h(jj,:),x(jj,:));   %% x is your signal

end

good luck


Understanding anti aliasing filter I'm working on a project for the university about the anti aliasing filter.

I have a small audio file, which has a sampling rate of 22KHz. Now when I downsample it 8 times, the file sounds weird. So far so good.

Now I should design an anti aliasing filter to avoid that. But how do I choose the cut off frequency based on the spectrum of the 8-time downsampled signal?

Anyhelp would be appreciated.

From your statement, i understand a discrete-time downsampling of a sequence $x[n]$ already converted to digital, without aliasing, previously.

So to prevent any aliasing on the downsampled sequence $y[n]$, you should first apply an anti-aliasing lowpass filter to the input sequence $x[n]$. The cutoff frequency of the lowpass filter is given by the downsampling ratio $D$ and is:

$$ \omega_c = \frac{ \pi }{D }$$ radians per sample. The gain of the lowpass filter is one.

The following Matlab\Octave code designs a linear phase, type-I, FIR lowpass filter of odd length $L = 2K+1$ with  a group delay of $K$ samples :

D = 8;     % downsampling factor
K = 256;   % group delay of FIR linear phase filter.
h = fir1(2*K, 1/D); % FIR, LPF, Type-I, L = 2K+1, wc = pi/D


P2 symbol constellation in DVB-T2 DVB-T2, 32K SISO, not extended.

I've acquired P2 symbol, PRBS bits have matched with bits recovered from pilots. Then equalization was applied to all subcarriers. The constellation I've got looks exactly like rotated 16-QAM overlapped with non-rotated 64-QAM.
On the picture: the red dots(-1,0) and (1,0) are for pilots, the blue dots for everything else, except tone reservation.

DVB-T2 Implementation guidelines (ETSI TS 102 831 V1.2.1 (2012-08)) clearly states:

  Constellation rotation shall only be used for the common PLPs and the data PLPs and never for the cells of the L1
  signalling.

What do I miss?
Also, can't find if there any subcarriers (cells) inverleaving applied to P2 symbol, can anyone explain this?

Ok, question solved. I found out that indeed L1 pre/post constellations are not rotated.

The constellation above comprises:

L1 pre (BPSK)
L1 post (64-QAM)
Data PLP from P2 symbol (rotated 16-QAM)

Also P2 cells are interleaved (ETSI EN 302 755 @ 8.5).


Choosing minimum sampling frequency of a filtered signal I have a time series of acceleration sampled at 1000Hz to capture random impacts on a structure. The ringing of the structure (not of interest) produces noise which a 50Hz filter seems to eliminate while preserving the impulses' amplitude and width. I'd like to downsample the filtered data to reduce data volume. What is the minimum sample rate I could choose? Based on this article, I should choose a sample frequency higher than 2*50Hz (Nyquist) since 50Hz is just the filter cutoff and there is still energy at higher frequencies that I don't want to alias. Is 10*50Hz = 500Hz reasonable?

As you say, you should select a sampling frequency at least at twice the largest frequency you wish to keep.

What you'll do is called downsampling: going from a higher to a lower sampling frequency. This process usually starts by filtering, so there is no energy above the new Nyquist frequency, and any aliasing is negligible.

Many numerical programs, such as Matlab, include the filtering step in commands such as downsample and resample. In any case, you should always look at the spectrum before and after downsampling, and make sure the process meets your requirements.


Best "Advanced DSP" book Let's suppose I have finished studying the first college level DSP textbook and I'm looking for a book, that is the next step beyond:
"Digital Signal Processing: Alan V. Oppenheim"

What's the best next book for "Advanced DSP" that would be best to study if on a self-study course?

(self-study meaning, it has sufficient examples and enough technical depth beyond an introduction textbook.)

Assuming you have completed Oppenheim's Discrete-Time Signal Processing book then the next (advanced) step could be anyone of the followings: (assumes a graduate level Linear System Theory background)

Statistical Digital Signal Processing, Monson Hayes
Multiresolution Signal Decomposition: Transforms, Subbands,Wavelets, A.Akansu
Adaptive Filter Theory, Simon Haykin
Multirate Digital Signal Processing, Rabiner
Discrete Random Signals and Statistical Signal Processing, C.Therrien
Fundamentals of Statistical Signal Processing : Estimation Theory, S.M.Kay
Multirate Systems and Filter Banks, P.P.Vaidyanathan
Wavelets and Filter Banks, Gilbert Strang
Multidimensional Digital Signal Processing, Dudgeon

To name just a few...


Analysing 2500 frequencies using FFT with an input vector of 2048 samples? I am currently reading the paper A Highly Robust Audio Fingerprinting System and on page 4 one can read about the technical parameters they use: Sampling rate of 5000 Hz, frames of 2048 samples as input to the FFT and with these settings they analyse frequencies from 300 Hz to 2000 Hz. Here are the respective quotes:

  Since  the  algorithm  only  takes  into  account  frequencies
  below 2kHz the received audio is first down sampled to a mono
  audio stream with a sampling rate of 5kHz.

.

  The most computationally demanding operation is
  the Fourier transform of every audio frame. In the down sampled
  audio signal a frame has a length of 2048 samples.

.

  These bands lie
  in  the  range  from  300Hz  to  2000Hz

What I don't understand: How can one analyse frequencies up to 2000 Hz using 2048 samples as input to the FFT? As far as I know the FFT output is mirrored, therefore with 2048 samples as input I get an output vector of length 2048, however what I can effectively really use are the first 1024 values (2048/2 due to mirroring). So I would have only 1024 values for 2500 different frequencies? (sample rate is 5000 Hz and according to Nyquist maximum possible frequency therefore 2500 Hz)

Thanks for any hint!

The first 1024 output bins of the FFT represent frequencies from 0 Hz to Fs / 2 = 2.5 kHz. So each bin is 2500 / 1024 = 2.44 Hz wide. In other words you have a resolution of 2.44 Hz in the frequency domain.


Where can I read up on zero state and zero input response for Discrete Time LTIs? There seems to be something I'm just not understanding about this topic.
Every question I work out that asks for zero state or input I get wrong.

Hence I think I am just not understanding the definition properly, could someone guide me to a preferably beginner friendly text that covers this topic in detail?

One of the Examples:

In 6b) I honestly have no idea how to tackle the question, but here's what I did:

I tried convoluting h[n] with u[n]. My reasoning is, that this way I can obtain y[n].
Then through some magic, I can get the zero state.

Here's how I derived y[n]:

Yes it's a little confusing at first.  And it only deepens if you also take an LTI system perspective into account, which imposes its own constraints into the corresponding difference equation structure.

Consider the following difference equation structure for a causal system:

$$\sum_{k=0}^{N} a_k y[n-k] = \sum_{k=0}^{M} b_k x[n-k]$$

with auxiliary conditions of $y[n]=0$ for $n<0$; known as initial rest.

The Zero State response of this system is equivalent to computing the output $y[n]$ for a given input $x[n]$ with initial rest assumption; i.e., setting all past values of $y[n]$ (aka the states of the filter) to zero before the signal is applied. This output $y[n]$ will equivalently be obtained by the well known convolution sum

$$y[n] = x[n] \star h[n] = \sum_{k} x[k]h[n-k]$$

where $h[n]$ is the impulse response of the associated LTI system described by the above difference equation under initial rest conditions.

The Zero Input response of an LTI system must be identically zero; which is evident by the definition of an LTI system : it must posses initial rest hence zero initial conditions and if you also set $x[n]=0$ identically, then the output must also be zero identically. So an LTI system can not have zero input response.

What happens in most of those questions is that, there is an LCCDE with arbitrary initial conditions; hence it won't strictly represent an LTI system, but it's output can be superposed into two components: one coming from input alone (assuming zero initial conditions) denoted as zero-state response, and the other coming from non-zero initial states alone, assuming no input but a response due only to the non-zero states of the recursive system.


Name of peak finding algorithm I've been using this implementation of a peak finding algorithm and it works pretty well. I'm wondering if this particular algorithm has a name?

Python code:

def indexes(y, thres=0.3, min_dist=1):
    if isinstance(y, np.ndarray) and np.issubdtype(y.dtype, np.unsignedinteger):
        raise ValueError("y must be signed")

    thres = thres * (np.max(y) - np.min(y)) + np.min(y)
    min_dist = int(min_dist)

    # compute first order difference
    dy = np.diff(y)

    # propagate left and right values successively to fill all plateau pixels (0-value)
    zeros,=np.where(dy == 0)

    # check if the singal is totally flat
    if len(zeros) == len(y) - 1:
        return np.array([])

    while len(zeros):
        # add pixels 2 by 2 to propagate left and right value onto the zero-value pixel
        zerosr = np.hstack([dy[1:], 0.])
        zerosl = np.hstack([0., dy[:-1]])

        # replace 0 with right value if non zero
        dy[zeros]=zerosr[zeros]
        zeros,=np.where(dy == 0)

        # replace 0 with left value if non zero
        dy[zeros]=zerosl[zeros]
        zeros,=np.where(dy == 0)

    # find the peaks by using the first order difference
    peaks = np.where((np.hstack([dy, 0.]) < 0.)
                     & (np.hstack([0., dy]) > 0.)
                     & (y > thres))[0]

    # handle multiple peaks, respecting the minimum distance
    if peaks.size > 1 and min_dist > 1:
        highest = peaks[np.argsort(y[peaks])][::-1]
        rem = np.ones(y.size, dtype=bool)
        rem[peaks] = False

        for peak in highest:
            if not rem[peak]:
                sl = slice(max(0, peak - min_dist), peak + min_dist + 1)
                rem[sl] = True
                rem[peak] = False

        peaks = np.arange(y.size)[~rem]

    return peaks

The full python package can be found here

The earliest use of the word "peak" in the IEEE library is from 1937

H. Roder, "Effects of Tuned Circuits upon a Frequency Modulated Signal," in Proceedings of the Institute of Radio Engineers, vol. 25, no. 12, pp. 1617-1647, Dec. 1937.
doi: 10.1109/JRPROC.1937.228819
Abstract: Prior investigations indicated that the frequency modulated receiver would always respond to the signal having the largest amplitude. Thus, selective circuits would be required to pick out a desired signal existing simultaneously with a number of other signals. The first item considered in this paper is that the signal carrier is tuned to the steep side of the resonance curve. It is found that in this case conversion from frequency modulation into amplitude modulation can be effected. It is required that the amplitude and phase characteristic of the circuit be linear with respect to frequency over the whole frequency interval occupied by the modulated signal. In order to derive a faithful audio signal reproduction from the complete detection process, a phase modulated signal must be transformed into frequency modulation first. Next, the case is considered in which the signal carrier is tuned to the peak of the resonance curve. The single tuned circuit is taken up first and two methods of solution are presented. It is found that with certain values of the modulation index, the tuned circuit may cause very serious nonlinear distortion of the output. For large values of the modulation index the resulting variation in amplitude and frequency can be determined statically, while for very small values of the modulation index the effect of the tuned circuit is exactly analogous to that encountered with amplitude modulation.
keywords: {Amplitude modulation;Chirp modulation;Frequency modulation;Nonlinear distortion;Phase detection;Phase frequency detector;Phase modulation;RLC circuits;Resonance;Signal processing},
URL: http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=1686321&isnumber=35536

As one of the commentators stated,  it's a very basic operation but also one that is not trivial.  The definition that I've used for myself is that a peak is a cell or bin that is higher than it's adjacent neighbors. It sometimes is the maximum likelihood estimate of discrete valued signal parameter.  This really ages me, but on a 300 pound, Numerix, Mars 432 array processor, searching through an array and testing if the left and right values of an element of a vector for "peakness" was a slow operation and it was faster to do two off set vector subtractions, clip and do a vector "and".


ZP-OFDM vs. CP-OFDM channel equalization Now, I am working on OFDM equalization, according to what I understood, Adding CP into the OFDM block change the convolution from linear to circular. That circular convolution enables using point-wise equalizer or one tap equalizer such that ZF or MMSE equalizer.
ZP-OFDM is a different scheme, where the CP is replaced by zeros. As I understood, ZP-OFDM should achieve the same performance of CP but under more complex equalization.
So, What means by complex equalization?? Is it time-domain equalization ? Could someone explain for me which equalizer we can adopt with ZP-OFDM in order to achieve similar/near performance of CP-OFDM ?

Hint: you could convert zero padded ofdm into cyclic ofdm at the receiver, and use the same equalization. Though there will be more noise on this scheme.
I had asked a puzzle some time back on this topic:
Modelling zero padded OFDM as circular convolution between channel and FFT frame


Recommended C/C++ General Purpose DSP Libraries I am looking to transition some of my MATLAB algorithms to C/C++ for deployment on target (PowerPC and/or ARM).  Mostly 1D FIR/IIR filter in this instance.

In other domains, there are well-regarded general libraries (both speed and quality) with liberal licencing. For example, Linear Algebra has LAPACK and Eigen (amongst others), and computer vision has OpenCV.

There are well-regarded libraries for specific algorithms, such as FFTW, but what about a collection of general algorithms for 1D signals?

I'm not after DSP-specific libraries, such as for PIC32, CMSIS (ARM Cortex M4), DSPLIB (C64x) or Renesas, but nonetheless has similar functionality.

Maybe a bit late, but since others might land here like me.

The following are good signal processing libraries/frameworks:
  * https://www.gnuradio.org/ - lots of classic signal processing + cascading options
  * https://root.cern.ch/ - advanced statistical signal processing
Both provide generic facade, fall-back implementations as well as performance optimised sub-libraries (ie. fftw, atlas, LAPACK, gsl, VOLK, ...).

Together with OpenCV, these probably cover 90++% of use-cases.


How is digital signal processing related to Computer Vision? Keeping in reference this image,

I want to build a concrete understanding of interconnection between Digital Signal Processing and Computer Vision. What would be the applications of this co-relation?

My attempt: You can treat an image as a discretization of the perspective projection of a continous 3D region, given a near and far plane. Using this discretization, filters used in signal processing can be used to modify the signal (a 2D array) for various purposes. These purposes include things like edge detection, image denoising, object detection, instance/semantic/general segmentation, and various other tasks. Notably 2-D wavelet transforms borne out of signal processing are used to transform image signals to another basis using function spaces spanned by a finite basis of orthogonal functions and perform filtering in this space, after which you invert the transform to achieve a modified image.


Proof regarding the periodicity of a continuous-time sinusoid after sampling Question
A continuous-time sinusoid $x_a(t)$ with fundamental period $T_p = \frac{1}{F_0}$ is sampled at a rate $F_s = \frac 1 T$ to produce a discrete-time sinusoid $x(n) = x_a(nT)$.

Show that $x(n)$ is periodic if $\frac{T}{T_p} = \frac K N$ (that is, $\frac{T}{T_p}$ is a rational number)
If $x(n)$ is periodic, what is the fundamental period $T_p$ in seconds?

I feel there is some ambiguity in the question itself. I assume $N$ here refers to the period of the sampled discrete-time signal.

I know how to prove that a discrete-time signal is periodic only if its relative frequency is a rational number. But the concept of periodicity in both continuous-time as well as sampled signal confused me.

So here is what I have tried:

Consider a continuous time sinusoid, $$x_a(t) = \cos(2 \pi Ft) \tag{1}$$ For this signal to be continous with a fundamental period $T_p$, it must follow that $$F=\frac{k_1}{T_p} \tag{2}$$ for some $k_1$ which is an integer.

Now that this signal is sampled with samppling rate $F_s$, the resulting discrete-time signal will be:
$$x(n) = \cos(2\pi f n)\tag{3}$$ where, $$f = \frac{F}{F_s} \tag{4}$$ This sampled signal will be periodic with period $N$ iff
$$f = \frac{k_2}{N} \tag{5}$$ for some $k_2$ which is also an integer.

Using equations (2), (4) and (5), it follows
$$\frac{T}{T_p} = \frac{k_2}{k_1 N} \tag{6}$$ If I redefine $\frac{k_2}{k_1}$ to be $K$, then equation (6) becomes
$$\frac{T}{T_p} = \frac{K}{N} \tag{7}$$ which proves the first part.

Now that the fundamental period $T_p$ has been asked, solving for $T_p$ in equation (6) gives
$$T_p = \frac{k_1 T N}{k_2} \tag{8}$$

Since both discrete time period and continuous time periods are in the same place, this has made me confused whether my approach to the solution is correct or not. Could you please verify that the approach is correct? Or if anything is wrong, could you point that out to me? Thanks.

You might be making this problem harder than it is by thinking in terms of frequency rather than periods.  Your sinusoid is
$$x(t) = \cos\left(2\pi \frac{t}{T_p}\right)$$ and its samples are
$$x[n] = x(nT) = \cos\left(2\pi \frac{nT}{T_p}\right).$$
The samples form a discrete-time periodic sequence of period $N$ if
\begin{align}
x[n] &= x[n+N]~\forall\, n\\
&{\Large \Downarrow}\\
\cos\left(2\pi \frac{nT}{T_p}\right) &= \cos\left(2\pi \frac{(n+N)T}{T_p}\right)\\
&{\Large \Downarrow}\\
\cos\left(2\pi \frac{nT}{T_p}\right)&= \cos\left(2\pi \frac{nT}{T_p} + 2\pi \frac{NT}{T_p}\right)
\end{align}
where the last equality will hold only if $\displaystyle\frac{NT}{T_p}$ equals some integer $k$. This is equivalent to the OP's desired requirement that $$\frac{T}{T_p} = \frac kN$$ for some integers $k$ and $N$.

So, what is the fundamental period of $x[n]$? That is,
what is the smallest integer $N_0$ such that $\displaystyle\frac{N_0T}{T_p}$ is an integer?  Well, suppose that the ratio $\displaystyle\frac{T}{T_p}$ is the rational number $\frac{k_1}{k_2}$ where $k_1$ and $k_2$ are relatively prime integers. Then, $N_0 = k_2$ and $\displaystyle\frac{N_0T}{T_p} = k_1$.  Thus, considering $x[n]$ as a periodic discrete-time sequence of samples spaced $T$ seconds apart  -- think of it as the impulse train $\sum_n x(nT)\delta(t-nT)$ -- gives a (continuous-time) period of $k_2T$ seconds for the sequence/impulse train. Note that one period of the discrete-time sequence spans $k_1$ periods of the sinusoid of period $T_p$.


What do we mean by spectrum? Consider the context of image processing and computer vision, and, in particular, discrete Fourier transform.
For example, in the sentence

In the Discrete Time Fourier Transform the forward Fourier Transform correspond to a discrete function of a sequence $x[k]$, however the inverse transform still remains continuous:
$$x(t) = \frac{1}{2\pi} \int_{-\pi}^{\pi} \hat{x}(\omega) e^{j \omega t} d \omega, \; \forall t \in \mathbb{R}$$
the inverse Fourier Transform is typically limited to integration on $[−\pi, \pi]$ (with $T = 1)$ as frequencies outside of the interval just correspond to replicas of the original spectrum produced by the sampling procedure.

or in the sentence

To make the inverse transform treatable by modern digital computer, we need to discretize the spectrum of the signal as well

The word "spectrum" is used. I have seen it being used in several other places. However, I still do not get its meaning.
What is the spectrum? I think this is related to the concepts of time, frequency, spatial and spectral domains, but how exactly?
I have a little understanding of Fourier transform, but I haven't yet fully grasped the concept of discrete Fourier transform. Furthermore, my knowledge of signal processing, image processing and computer vision is very limited.

The "spectrum" refers to the frequency content of the signal (both phase and amplitude/power). For a time series or 1 dimensional signal, it basically represents the presence or non presence of the different possible gradients in the time domain signal.

For an image the spectrum simply represents how quickly/slowly the pixels change in contrast/color/values in both the spatial dimensions.


Symbol time offset simulation I ran a number of simulations for symbol timing synchronization in the past. I used GNU Radio mainly. Now, I'm trying to do the same thing in C++. To achieve that, I need a way to introduce symbol offset to the symbols. So suppose I have a complex array of $M$ samples acquired at a sampling rate of $N$ samples per symbol. How do I introduce a symbol offset of say $\mu = 100 \ \rm ppm$ of the symbol time? Pseudocode, C++, or MATLAB script will suffice.

If you postulate that receiever's clock is perfect, then you want to make the transmitter send symbols every $T_s \pm \varepsilon$ seconds, where $T_s$ is the symbol period according to the receiver.
This is easily achieved by using a very high sampling rate in the transmitter. Let's assume $T_s=1$ and you need a deviation of $\pm 0.01$. This deviation corresponds to one sample, if the sampling interval is $0.01$. So, use a sampling rate of $1/.01 = 100$ samples per time unit, and have the transmitter send a new symbol every 99 samples (for a faster clock) or every 101 samples (for a slower clock). Then, have the receiver sample every 100 samples, and you're done.
By making the sampling rate fast enough, you can easily accomodate any delay you need. Using this method, it is also easy to simulate a time-varying clock drift in the transmitter and/or receiver.


Need a 2 parameter ranking system for image compression I'm working on a project in which specific bit planes of an image are extracted and compressed. This process is repeated for every combination of bit planes possible.For an image, I have calculated the compression ratio and the Structural Similarity Index(SSIM).

I need a method to allot ranks when compressed with different combinations of bit planes.

My progress:
I tried to give equal weight to the compression ratio and the SSIM value.

score = 0.5*ssim + 0.5*compression_ratio;

This equation gave me a value between 0 and 1. The problem with this equation is it gives an equal score for multiple combinations and few combinations are extreme opposites on the scale. For example, i got equal scores when I compressed an image using the {1,2,3,5,7,8} bit planes(High SSIM value) and only the 3rd bit plane(High compression ratio).

I need help in balancing the parameters and getting an appropriate score. The score doesn't need to give equal weights like I have given.

Thanks in advance

I am not sure if this is going to end up having 2 parameters specifically, but there are a number of things to note:

score depends on two parameters right now, but it is clearly unable to resolve the differences between different schemes. You can add a third parameter which captures the bit planes you are using. This could be either via:

Normalising the number of bit planes used (e.g. $\frac{3}{8} = 0.375$)
but this might still be not enough as a simple reshuffling of the bitplanes can show;
Indexing your score by the number that results from the bitplanes used. So, if you used bitplanes 0,3,5, that would be 2^0 + 2^3 + 2^5, or simply a tuple of $(0,3,5)$ so that implicitly you have the number of bitplanes and which ones were used. This would give you a lookup of (bitplanes_used, score).

Now, the danger you are running here is formulating some $score = f(\Theta)$, where $f$ is a function whose parameter set $\Theta$ keeps growing, more than the 2 you specify. What if we include the ratio of high frequencies to low frequencies in there too? What if we add an entropy related term?...and so on. Ultimately, you might end up with a very large model that is still not informative.

A common thing to do in this case is to penalise large models. This is an extra factor that increases with the number of parameters used. For example, in your case, you could say something like "If two bitplane combinations give me the same score, then I will select the combination that uses fewer bitplanes than the other"

But, it is possible that after all this, your lookups still contain equal scores. In that case, you can collect a subset of the values (that satisfy some criterion in your application) and then randomly select one of them, since all of these values represent a choice of equal score.

Hope this helps


FFT function isn't working on imported signals greater than a second I have an fft function in matlab/octave that works great on an imported mono wav/audio signal that is 1 second long, but when the signal imported is more than 1 second the frequency values become incorrect.  To check this I created a 2 second signal using audacity the first half of the 2 second signal (1 second of it) was 500hz at 0.3 amplitude and the second half of the 2 second signal (1 second of it) was 200hz at 0.8 amplitude.  I also checked it with audacity's analyse spectrum option.

When I use the function on this 2 second signal it shows that the max frequency is at 999.95hz when it should show the max frequency is at 200hz.  Can anyone help me get this function to work with signals greater than just 1 second.

I've included the function along with the code that calls the function.

The FUNCTION:

function [freq,amp,ampinv,phase,phase_radtodeg,phase_degtorad,phase_adjindeg,phase_adjinrad,sigfft,sigifft,sigphase]=rtfftphase(vp_sig_orig,Fs,phase_goal)

    vp_sig_orig = vp_sig_orig - mean(vp_sig_orig); %remove the mean of your signal due to dc offset

    vp_sig_orig=vp_sig_orig';
    vp_sig_len=length(vp_sig_orig); %get sample rate from vp fs_rate needs to be an even number?

    % Use next highest power of 2 greater than or equal to length(x) to calculate FFT.
    nfft= 2^(nextpow2(length(vp_sig_orig)));

    % Take fft, padding with zeros so that length(fftx) is equal to nfft
    fftx = fft(vp_sig_orig,nfft);
    sigfft= fft(vp_sig_orig);
    sigifft=ifft(sigfft);
    sigphase = unwrap(angle(sigfft')); %get phase of orginal signal

    % Calculate the number of unique points
    NumUniquePts = ceil((nfft+1)/2);

    % FFT is symmetric, throw away second half
    fftx = fftx(1:NumUniquePts);

    % Take the magnitude of fft of x and scale the fft so that it is not a function of the length of x
    mx = abs(fftx)/length(vp_sig_orig); %replaced for testing from stackexchage

% Since we dropped half the FFT, we multiply mx by 2 to keep the same energy.
    % The DC component and Nyquist component, if it exists, are unique and should not be multiplied by 2.
    if rem(nfft, 2) % odd nfft excludes Nyquist point
    mx(2:end) = mx(2:end)*2;
    else
    mx(2:end -1) = mx(2:end -1)*2;
    end

    %yamp=(mx(1,:)/max(abs(mx(1,:)))*1); % keep at 1, dont use in function use in script file directly
    amp=mx;
    ampinv=abs(amp-max(amp));

    % This is an evenly spaced frequency vector with NumUniquePts points.
    freq_vect = (0:NumUniquePts-1)*vp_sig_len/nfft;
    freq=(freq_vect'); %take of round if you want exact numbers

    %get phase of new signal
    phase = unwrap(angle(fftx)); %get phase of orginal signal

    %phase stuff phase adj stuff not fully done yet I don't know the correct way to use it leave at 0 for now
    phase_radtodeg=rtrad2deg(phase); %also check sigphase to see if it's the same.
    phase_degtorad=rtdeg2rad(phase_radtodeg); %for checking
    phase_adjindeg=phase_goal-phase_radtodeg; %used to create new phase to get to goal phase
    phase_adjinrad=rtdeg2rad(phase_adjindeg);

Example: that calls the function

imported_sig_L=imported_sig(:,1)';
[xfreq_orig,yamp_orig,yamp_inv,phase_orig,phase_in_deg,phase_in_rad,phase_adjindeg,phase_adjinrad,sigfft,sigifft,sigphase]=rtfftphase(imported_sig_L,fs_rate,0);

array1=[xfreq_orig,yamp_orig,yamp_inv,phase_orig,mod(phase_in_deg,360),phase_in_rad,phase_adjindeg,phase_adjinrad]; %have degs only go up to 360

%Sort Array Section

array1_sort_amp_norm=sortrows(array1,-2); %sort by amplitude
array1_sort_amp_inv=sortrows(array1,-3); %sort by amplitude inverse

%get max freq checking to make sure freq isn't zero if so go to next row
if (array1_sort_amp_norm(1,1)>0); % To get max amplitude
    maxfreq_normA=array1_sort_amp_norm(1,1)
elseif   (array1_sort_amp_norm(1,1)<=0); %
    maxfreq_normA=array1_sort_amp_norm(2,1)
end

if (array1_sort_amp_inv(1,1)>0); % To get inverted max amplitude
    maxfreq_invA=array1_sort_amp_inv(1,1)
elseif   (array1_sort_amp_inv(1,1)<=0); %
    maxfreq_invA=array1_sort_amp_inv(2,1)
end

Thanks

Here a short example of how to deal with FFT in octave/matlab:

function freq(npts, len, basefreq)
    x = [0 : npts - 1] / npts * len;
    y = sin(basefreq * x * 2 * pi);
    F = abs(fft(y));
    half = npts / 2;
    F = F(1 : half);
    plot([0 : half - 1] / len, F);
    printf("Base frequency = %g\n", (find(F == max(F)) - 1) / len);
endfunction

try this function with different parameters.

Here npts is size of data vector; len - time length of data (in seconds); basefreq - base frequency for sine function.
First we fit a time axe (x), after that we calculate sinusoid with given frequency, get modulus of its FFT, trim it to half (BTW, another variant is not to trim, but to move frequency zero to middle of fftshift'ed data). After all we plot result of FFT and print out founded base frequency.


Analyzing the quality of a music track I have a library of music tracks that I use to DJ with. Its currently about 3000 tracks that I have gather over the years. Some of it consists of low quality rips that I want to get rid of.

Currently I am writing a script that will look at the time/size = compression rate. Anything that is below 300kb/s I want to delete.

I am wondering if that is a good heuristic to use, or could I possibly do something fancier. Like looking at a spectral analysis to determine what is low vs high quality.

Currently I am writing a script that will look at the time/size = compression rate. Anything that is below 300kb/s I want to delete.

That's a very bad idea.

Just because something compresses well, because it fits the signal model of the compressor well, doesn't mean it's bad quality. Usually, quite the contrary.

You give an example in your comments yourself:

  Im talking about tracks that have lot of clipping

Clipping is a nonlinearity, introducing a multitude of tones into the spectrum that shouldn't be there. Since psychoacoustic models and thus lossy compressors usually need to quantize the spectrum in some way, having clipping in your audio means harder to compress audio meaning larger file size.

Generally, if something is hard to compress, a variable-rate codec will increase the file size. If something is easy to compress without much loss, then the file size will stay down. I'm pretty sure you can record a pretty perfect tuning fork tone in high quality and compress it very nicely.

Also, unless we're talking early MP3 encoders and such, modern codecs above ca. 160 kb/s are really hard for a human to tell from the uncompressed audio. That's the whole idea of having a good lossy compressor: Just lose the information that's irrelevant.

Note that I'm not saying you can't generally say that for example an MP3 with 64 kb/s will sound bad in almost all cases. It's just that you don't need any intelligence to do that, but can often directly read that from the file's metadata.

  Like looking at a spectral analysis to determine what is low vs high quality.

Nope, you're looking at the reconstructed signal after decompression. Unless you have extensive knowledge of the spectrum of your original audio, there's little information about audible losses in that. Again, both musicians, sound engineers and audio codecs strive to do the same: produce sound that sounds good to the human ear.

Also, remember that clipping (or soft clipping, and other similar nonlinearities), as the example you chose, is actually a pretty commonly used tool in production of music these days. Who are you to tell someone that the recording of their neighbor's cat walking on a keyboard is quality-wise better than let's say a Kraftwerk song that uses clipping intentionally? How is a crisp-sounding KPop-girlband song with strongly emphasized treble quality-wise superior to the muddy sound of the 1968 original release of Jimi Hendrix' All along the watchtower? There's whole musical genres that depend on sound mechanics that would signify a lack of recording quality in others. Even within the same genre, quality can't simply be deducted from a rigid description of the signal – compare Shostakovich's The Bolt to other ballet suites of that era.

"Good audio quality" in the end is an extremely subjective matter, and I'm afraid without a large-scale database of manually tagged "good" and "bad" examples according to your perception, you won't be able to implement any classifier, be it through neural nets or through more classic statistical approaches.


Is $x[n]=(-1)^{n^2}$ periodic? Is $x[n]=(-1)^{n^2}$ periodic? The answer said no, but when I draw it on a graph, it seems to be periodic, with fundamental period equal to $2$.

You're right, the given $x[n]$ is clearly periodic. You can show this by simply checking if

$$x[n]\stackrel{?}{=}x[n+N]\tag{1}$$

is satisfied for some positive integer $N$.

For the given $x[n]$ you get

$$\begin{align}(-1)^{n^2}&\stackrel{?}{=}(-1)^{(n+N)^2}\\&=(-1)^{n^2}(-1)^{N(N+2n)}\tag{2}\end{align}$$

From $(2)$ it follows that for $(1)$ (and $(2)$) to hold we require

$$(-1)^{N(N+2n)}=1\tag{3}$$

i.e., $N(N+2n)$ must be even (for arbitrary integer $n$). This is certainly the case for any even $N$, and the smallest positive even $N$ is $2$, hence the (fundamental) period of the given $x[n]$ is indeed $2$.


Wavelets: How does the admissibility condition imply this $\Psi(0)=0$ I'm reading through some material on wavelet and I'm finding it difficult to see why the following is true.

Let $\Psi(\omega)$ is the Fourier transform of the mother wavelet $\psi(t)$

$$\int_{-\infty}^{\infty}\frac{|\Psi(\omega)|^2}{|\omega|}d\omega<\infty$$

then many paper on wavelet would immediate jump to, is obvious that:

$\Psi(0)=0$

But wouldn't the above yield a 0/0 condition?

I've also tried relaxing the limits and take a derivative on both side, but this would only imply that the mother wavelet is zero at some frequency $\omega$ not $\omega$ = 0. How do you rigorously prove the above relation?

It's obvious that Fourier transform of the mother wavelet is zero at the origin, otherwise the integral, whose argument is nonnegative, would be unbounded. Just because both the numerator and the denominator are zero, it does not mean it does not have a limit. So we know that $\Psi(0)=0$, but this is merely necessary, not sufficient for the integral to be bounded.

Consider the function $f: \omega \to |\sin(\omega)|^2/|\omega|$. How does this function's integral from $(-T,T)$ behave as $T \to \infty$. I hope that gives you some ideas.


Output of discrete-time LTI system guaranteed to be same form as input? I know that in the continuous-time context, if I supply a complex exponential input to a Linear Time Invariant system, the output will be of the same form as the input - for example, if the input is $x(t)=Ae^{j \omega_0 t+\phi}$, the output will be of the form $y(t)=|H(j\omega_o)|\centerdot A \centerdot e^{j(H(s)\omega_0 t+\phi +\theta_H)}$ where $H(s)|_{s=j\omega_0} = |H(j\omega_0)|e^{j\theta_H}$ is the transfer function of the system evaluated at the frequency $j\omega_0$. My point is that the system does not modify the input frequency.

Does the same hold true in the discrete-time context? If I supply an input $x[n]=a^n$ to an LTI system, will the output also contain $a^n$? Or is it possible to supply $a^n$ and get the output $b^n +c^n$ where $a\ne b\ne c$ (for example, $x[n]=(1/6)^n$, $y[n]=B(1/2)^n + C(1/3)^n$)?

Edit: I've posted a related question that contains the issue that actually prompted this question...

The point is that a complex exponential function is an eigenfunction of an LTI system link. Therefore the frequency is not changing, i.e., the LTI system translates to a (complex-valued) gain $H(j\omega_0)$.

For discrete systems, the eigenfunctions are sampled exponential functions link, i.e.,
$x[n] = A\cdot e^{j\omega_n n}$, where $\omega_0=2\pi f/f_s$ is a normalized frequency. The output of a discrete LTI system is then $y[n]=H(e^{j\omega_n})\cdot A \cdot e^{j\omega_nn}$.

EDIT:
A discrete input signal x[n] can be written as a weighted sum of complex exponentials, i.e., $X(e^{j\omega}) = \sum_{n}x[n]e^{j\omega}$ link. This transformation is in fact a change of basis. Driving such a signal thru an LTI system will weight each exponential function with the complex-valued gain $H(e^{j\omega})$. Thus, no additional frequencies will be introduced.

Note (a bit out of scope, but good to realize):
When dealing with sampling a contuous-time signal and dealing with finite length signals, then you may observe additional frequencies due to aliasing and windowing.


Comparing the control performance of a given closed loop system I am having a question about comparing the control performance of a given closed loop. In this link there are some characteristic values that need to be calculated, but there is no description of what they mean or their interpretation for example: the root-mean squared control error and root-mean squared change in the manipulated variable (manipulating effort).

This is the formula for the root-mean squared control error :

$$
S_e = \sqrt{\frac{1}{M} \sum e^2(k)}
$$

Where the $e$ is the error

And root-mean squared change in the manipulated variable :

$$
S_u = \sqrt{\frac{1}{M} \sum \Delta u^2(k)}
$$

where

$$
\Delta u = u(k)-u(\infty)
$$

Could someone please explain their source or their interpretation ?

The error is a system state's distance from a desired state. If you want a dart to hit bullseye on the dartboard, and you miss, the distance to the bullseye is the error. The RMSE (root-mean-square error) is a way to summarize the error over time, expressing the error as one number. This makes it more easy to compare different cases. For a stationary/steady-state reference signal, such as a sinusoid, if you subtract the mean, it converges to the standard deviation of the error. The sensitivity plot for a closed loop system gives the RMSE directly for a given frequency. Think of a reference signal frequency, look it up on the Bode diagram for the sensitivity, and you know the scaling of the error at that frequency.

The root-mean-square change is a way to measure the power spent on controlling your system. Often you have limited resources available to drive a system to a given state. Think of running vs. walking. If you can slow things down, you reduce the power required. If you draw too much power you can exceed the limitations of a physical system. You might saturate a voltage amplifier or snap an axle... There is a classical trade-off between manipulation power and the smallness of the error.


How can detect the Left Ventricle (LV) Cavity into the MR dicom file? This is MR Modality Link

MR Modality

In this file where i find the LV Cavity in above link file. Also how to draw Tumor Marker from over it

I'm looking at this DICOM now and it is very hard to tell where I am in the brain (since there is only one slice). Anyhow, tumours tend to be bright irregularities. The thing I think is a tumour certainly isn't in the LV.

You can just draw an ROI over it using something like FSLVIEW if you are planning to do it manually, or you can have some fun and use something that looks for irregularities across hemispheres in the edges as detected by an algorithm.

http://www.medicalimagingandgraphics.com/article/S0895-6111(11)00079-6/abstract

Hope this helps...


How to design a delay and sum beam former with fractional delay less than one sample? I am supposed to design a delay and sum beam former for a short audio. But the sampling frequency is only 8kHz and after times the delay, the fractional delay is less then one sample, and after the "floor" function in Matlab, all the fractional delay integer will be 0 and the results are neither fractional delayed nor time aligned. Exponential weights cannot be used in this case for the reason no imaginary part should be in the filtered signal.

There is no reason you cannot use exponential weights. Those are used in the frequency domain, not in time domain, and won't produce imaginary part.

Using the Translation/ Time-Shifting Property of the Fourier transform, you can define the weights (per frequency bin), even for fractional delay. By the way, the problem of fractional delay is not only when it is less then 1 sample. Flooring all delays to be integers will produce sub-optimal results.

A different approach is to use interpolation in time-domain to obtain values between samples, e.g. using interp1. This approach will be much slower, if it matters.


Designing a lowpass filter to minimize aliasing in pre-decimated streaming audio I need to apply low-pass filter to PCM files. There are several methods such as FIR filters, IIR (butterworth-chebyshev..) filters but it seems to me applying Fast Fourier transform and eliminating higher frequencies is the closest way to an ideal filter.

What is the fastest and closest to ideal filtering method?

It is required for anti-aliasing before changing the sampling rate of sound.(Fcutoff=Fs/2) And will be applied every block of 1 sec. data. The main requirement is, after changing the sample rate, the new audio quality must be close as possible to original audio quality. (Not noisy.)

Thanks.

The best choice of filter depends on your specific application requirements. There are two basic choices: FIR and IIR. IIR will be much more efficient however, it will results in phase distortions. The phase distortions are completely inaudible (unless it's a bizarre outlier case) but clearly measurable. So it depends whether you can tolerate this our not.

In either case you need to decide how close you need to get to the new Nyquist frequency and how much aliasing noise you can tolerate. A typical example would be that you want the passband to extend to 90% of the new Nyquist frequency and that you would like your aliasing products to be below -80dB. Based on these specifications you can then design the appropriate filter. Other considerations include how much pass band ripple you can accept and if you have any constraints on maximum group delay and/or latency.

Here is an example: Let's say you want to downsample from 44.1 kHz to 32 kHz and the new Nyquist frequency is 16kHz. Going to 90% Nyquist (14400 Hz), with 0.1dB pass band ripple and 80 dB of attenuation at 16 kHz could be done with an elliptical filter of 9th order.

As nibot has pointed out, zeroing FFT bins is a poor choice for a low pass filter since the resulting low pass has very big side lobes and aliasing rejection will be quite poor. It would also require a proper implementation of an overlap-add or overlap-save algorithm to deal with a continuous signal.


In pulse code modulation what is pulse word represent? we study lot of things in PCM we do sampling quantization encoding etc but we never come across pulse so what is pulse word represent in PCM

Think of a DAC as a thing that is emitting a train of pulses, each new one hard up against the last one, and each one exactly as wide as a sampling interval (sampling interval = $T_s = \frac{1}{f_s}$).
The code in PCM is the height of each pulse.
And if the name seems odd now, just remember that it made perfect sense in the 1925 to 1948 time period when the basis for this stuff was being laid down.


Interpreting Cepstrum of plucked string in my attempts to synthesize a plucked electric guitar string (using supercollider's digital waveguide, which is based on the Karplus-Strong algorithm) I came to a dead end and decided to look that the Cepstrums of a sampled string and my synthesized string, to nail down the audible differences. They look like this (A-string - 110Hz)

My questions are:

The sampled signal has something in the right half, which the synthesized signal does not. I cannot see how to interpret this as "excitation". So what does it stand for?
Is it safe to say, that the spectrum of the synthesized spectrum is more "harmonic" than the one of the sampled signal? When I look at the spectra themselves it does appear this way. The sampled spectrum (left) looks more "noisy")

Could it be that the regularity (and the dull sound) of the synthesized signal is a consequence of the digital waverguide itself. The way I understand Karplus-Strong, after some time only the fundamental and its harmonics survive, which would indeed lead to a "harmonic" spectrum.
If that is the case, what could be the cause of the irregularties of a real plucked string. IOW: what does a real string have, which Karplus-Strong doesn't have?
As a side note: what does the sharp peak at the end of the Cepstrum stand for? Is it caused by the 1 second window?

IIRC, Karplus-Strong is usually a single linear harmonic waveguide.  Plucking a stringed instrument starts several things in vibration other that the string (sound board, air in cavity, etc.), many of them producing inharmonic overtones due to material stiffness, thickness and non-uniformity, some of which exchange energy and/of change in a non-linear fashion over time. The sharp peak high in the Cepstrum could indicate some kind of tiny but sudden frequency modulation(s) of the fundamental pitch.


Calculation number and meaning of DFT I'm posting this question to both check my own understand of DFTs (I'm relatively new to the subject), and to ask a few questions at the bottom regarding the DFT and its output. If anything isn't allowed then just let me know.

If we have a continuous signal in the time domain, $x(t) = cos(2\pi f t)$, we need to convert it to a discrete time signal.

If we sample this signal at a rate $f_s$,  which is the no. of samples taken per second, our equation becomes:
$$x[n] = cos(2\pi f n T_s)$$
$$\Rightarrow x[n] = cos(\frac{2\pi f n}{f_s})$$
where n is the sample number.

If we wanted to find the spectrum of this sampled signal, we could use the DTFT which would give us a continuous and periodic output. But because we're working with computers, the output must be finite, so instead we use the DTF which will give us a sampled version of the DTFT output. The DFT is given by:
$$X_k = \sum_{n=0}^{N-1} x_n . exp^{-j \frac{2\pi kn}{N}}$$
where N is the total number of samples, n is the sample of the input signal x[n], and k is the number of the frequency bin.

The frequency bin resolution is $\Delta f = \frac{f_s}{N}$, where N is the total no. of samples of x[n].

If we were to sample an input signal $x(t)$ with maximum frequency $f_{max}$, then by the Nyquist Sampling Theorem we want the sampling frequency $f_s = 2f_{max}$ to sample up to $f_{max}$ without aliasing.

I was wondering if my understanding of the DFT above is correct, and does this mean that, for each DFT bin, the DFT performs N summations and N multiplications to give the amplitude of the frequency at that DFT bin? And lastly, if I were to plot the amplitude vs frequency of the output of the DFT, would I only consider half of the spectrum (up to $\frac{f_s}{2}$) and then just double the amplitudes of the frequency bins below that?

Plain DFT performs  N summations and N multiplications to give the amplitude of the frequency at that DFT bin but usually, we use FTT which is an algorithm for reduction of the number of multiplications.
Since you project your signal on the complex exponent, your output is complex too. There are different ways to represent the output. You can decompose the spectrum to amplitude and phase (Bode plot) or just ignore the phase and plot only the amplitude.

For a real input signal, the spectrum should be symmetric around zero and you can present the right plane only but for a complex input signal, you have to represent both sides of the spectrum. Regarding doubling the amplitude, I am not sure that it is mandatory but it is a common practice.


Autocorrelation of a noisy linear map I am interested in calculating the autocorrelation function of a linear map with some noise (model given below) but am slightly confused in doing so.

At first, I did not realize there were two definitions of the autocorrelation function: one where you subtract the mean of the time series data (often used in statistics), and one where you do not.

Let $x(n) = \frac{1}{2}\left(a x(n) + b + \xi\right)$ where $\xi$ is some Gaussian noise with mean zero and standard deviation $\sigma$. One can analytically calculate the autocorrelation function as follows:

$$
\langle x(n+1)x(n)\rangle = \langle \frac{1}{2} (a x(n) + b + \xi_n)x(n)\rangle\\
= \frac{a}{2}\sigma_x^2 + \frac{b}{2}\mu_x
$$

Where we have used that the noise is uncorrelated with the time series data by assumption.

Similarly,

$$
\langle x(n+2)x(n) \rangle = \frac{a}{2}\langle x(n+1)x(n) \rangle + \frac{b}{2}\langle x(n) \rangle \\
= \left(\frac{a}{2}\right)^2 \sigma_x^2 + \frac{a}{2}\frac{b}{2}\mu_x + \frac{b}{2}\mu_x
$$

And then generalizing,

$$
\langle x(n+\tau)x(n) \rangle = \left(\frac{a}{2}\right)^\tau\sigma_x^2 + \frac{b}{2}\mu_x\sum_{l=0}^{\tau-1}\left(\frac{a}{2}\right)^l\\
= \left(\frac{a}{2}\right)^\tau\sigma_x^2 + \frac{b}{2}\mu_x\frac{1 - \left(\frac{a}{2}\right)^\tau}{1-\frac{a}{2}}
$$

Where I've summed the geometric series to arrive at the final result.

Autocorrelations are computed using NumPy's fftconvolve function as acf = fftconvolve(x, x[::-1]) (i.e., convolve the time series data with itself reversed, for anyone not familiar with Python).

When I subtract the mean, I get an exponentially decaying autocorrelation function, as depicted in the plot below. This makes sense, given the first term in the analytically derived expression above. Clearly $\mu_x = 0$ after we have subtracted off the mean.

However, originally, I did not realize I was supposed to subtract the mean, and I found a linearly decaying plot as depicted below. By increasing the duration of the simulation, the line extends seemingly to infinity. Is there a simple explanation for this? It does not fall out of the analytical calculation, so I assume there must be something going on because it is a discrete time series? I am new to this, but I assume there is a simple explanation for this phenomenon.


Why can't DFT be used when samples are not equally spaced in time? I found the following comment here

  The DTFT can be used when the samples are not equally spaced in time, the DFT cannot

My initial thought was that this had to do with periodicity of the basis complex exponential in DFT, but then this issue is unrelated to the uniformity of sampling itself.

What is the correct explanation here ?

You can actually use the DFT on non-unformly sampled data. It is still an orthogonal bijective map and has no inherent mathematical issues.

However, interpreting the meaning of the resulting spectrum is much less intuitive than for evenly spaced data.

The reason is that the basis functions are no longer eigenfunctions of the time translation operator but of an operator that shifts samples, i.e. each sample gets shifted by a possibly different time offset. This is not a very meaningful operation in general and the eigenvectors of it are also not terribly useful.

As a result the bins of the spectrum do not correspond to a single frequency, but rather a single rate of samples. It is possible to associate a single bin with a sinusoid of time-varying momentary frequency though.

This can still be useful, mostly if your sampling intervals are not random but intentionally spaced differently to follow some kind of underlying structure of your data.

So in total, saying that you can't use a DFT on non uniformly sampled data is not accurate. You can, but you have to know exactly what you're doing and be careful with drawing conclusions.

EDIT:

Let me add some remarks about the relationship between DFT and DTFT in this context. The commend of Dilip in the question you refer to has probably the following background. The discrete time Fourier transform can be written as the continuous time Fourier transform of a series Dirac delta distributions:

$$ \mathrm{DTFT}\{x[n]\} = \mathcal{F} \{ \sum_nx[n]\delta(t-t_n) \}$$

where the $t_n$ are the sampling times associated with the sampled values $x[n]$. This relationship naturally generalises to a nonuniform sampling grid if you just assume that the $t_n$ are not given by $nT$ for a sampling period $T$. I'm not sure I would call the result of this a DTFT anymore, but it's mathematically valid to do. because the properties of the continuous time Fourier transform make sure that the result is well behaved and that the basis functions are still orthogonal and represent a single frequency. This only works if your time axis is all real numbers, because only then the complex exponential functions expanded into are guaranteed to be both complete and orthogonal.

For a finite time axis where the DFT would be applicable, this method does not work. You can theoretically define a map that has all the properties of the DFT on a uniform grid when applied to a non-uniform grid by resampling the data on a uniform grid, applying the ordinary DFT, and considering these two steps as one "nonuniform DFT transform". In this case you will in general loose certain mathematical properties like unitary/orthogonality however.


Digital Image Processing Framework i'm looking for a free/opensource solution for digital image processing with reasonable performances processing quite large images (e.g. actual smartphones have resolutions of 8-10Mpx) and being ready for fast-prototyping/research.

Octave would be a good solution as it provides a lot of statistical and matrix manipulation tools. However octave is too slow when running nested loops and the provided image package is too slow when displaying images.
Octave + MIJ would be a better solution but actually seems unusable and unstable because of slowness and because of the different datatype handling needs.
OpenCV has excellent performances but it's not really handy for fast prototyping because i actually would need something like a scripting language (no compilation, online code/variable modification and so on...) and some tools (like octave has) for plotting functions, graphs, histograms and make statistical calculations.

Any advise?
Thank you very much

My first and second choices are Python and IDL, not necessarily in any order.

Python is popular with scientists, engineers and artists, well-known, on all platforms, free and open.  You need Numpy and PIL (Python Imaging Library).  Python is more and more commonly used in many places where Matlab used to rule.

IDL from IIT-VIS (I think it's a different name now) is similar to Matlab and the like, but nicer syntax for multidimensional arrays. Very popular with astronomers, geologists, and various others.  It has flaws.  IDL costs $$ but luckily there's a free version named Fawlty Language. It's not open source, but is free (last time I checked.) I don't have the link handy, but google it.   There's also GDL which is a Gnu attempt to create an IDL-compatible interpreter, but last time I tried it, it wasn't so good.

I don't like Octave because it's clumsy for three dimensional and higher arrays.  Images are only 2D, but I often need to process a stack of them, and also may have RGB or other types of channels.


AGAST vs FAST evaluation In ~2006 the FAST corner detector was introduced. One of it characteristics is that you have to train a decision tree from some test image in order to learn the types of corners present and relevant in the test image. One of the consequences is that the images in which you will be applying the detector have to be somewhat similar in relation to the pose (perspective), global illumination, etc. to the test image used in the training phase.

Feature descriptors such as BRIEF and ORB make use of FAST or a modified version of FAST in their reference papers.

Then, in ~2010 the AGAST corner detection was proposed to address this issue. According to the authors, you should not have to retrain the detector on a test image (but you may if you really want to or think this is really necessary).

Then feature descriptors such as BRISK and FREAK showed up, using a modified version of AGAST in their reference papers.

My question is: Since AGAST is an evolution of FAST aimed at solving some of its issues, should someone simply use AGAST from now on when building a detection application? (I am not speaking of doing pure research on the topic of detectors)

Is there some paper proving or at least showing convincingly that you should now use AGAST? I know I could simply test them both every time I build a detection application, but honestly, if there is a relatively good point I should use AGAST, it would be time saving.

I am not aware of a direct comparison, other than the agast paper:
http://www6.in.tum.de/Main/Publications/Mair2010c.pdf

However, Agast is also based on AST, which is developed by the authors of FAST -Edward Rosten (et al.). It is just a better training. So it might be safe to assume the improvements, based on the original Agast paper.

Moreover, here is a tracker I implemented based on Agast features:

https://www.youtube.com/watch?v=rZWu4LxmsME

This works fairly well in many cases and was able to convince me on the power of Agast (Doesn't work as good with FAST features).


Melodic or tonal pitch contour tracking? What algorithms or methods can be used for the tracking of audio pitch contours?  (Such as in musical vibrato and glissando, or for tonal vowels in languages such as Chinese?). Will a PLL work even for pitched tones that have a weak or missing fundamental?  Or is there a more suitable tracking method?

I think a PLL would be hard to use for pitch detection in singing or speech. For the reasons you mention (classic phase/frequency comparators output garbage when the waveforms are too different due to a missing f0), but also because you'll need another method anyway for tracking the large changes in pitch or the switch to voice/unvoiced that would put a PLL off the tracks.

A baseline system could be to extract a pitch on small frames with a classic f0 detector (autocorrelation, AMDF, YIN, or any auditory model), and then use dynamic programming (or from another point of view, the Viterbi algorithm on a HMM whose states are discretized pitch values) to recover a pitch track. A ready to use Matlab implementation for this (auditory model + viterbi) is available here or (fancy autocorrelation flavour + dynamic programming) here.

This will be discretized on both the time and pitch axes, but you could then interpolate.


How to calculate particles sizes for spheric particles with overlapping and superposition? (Example image included) I'm trying to automatically calculate particles sizes histogram for images like this one:

Common threshold and segmentation doesn't seem to work. What filters and algorithms do you suggest to be used?

Here is what I experimented with:

Use ELSD to generate elliptic contours. You could basically use any edge detector, but since in the following stages I will benefit from circle detectors, it is good to already have some geometrical edges. Here is what the output looks like:

Collect all circular (also non-circular) edges via collecting the channels into a single image. This can also be done via simple thresholding. Again, the output looks like:

I will denote this image to be $C$.

We will now compute the circular hough transform on this image. However, there are many spheres (which appear to be circular) in varying radii. Thus, we select a range $\mathbf{r}=\{9, 12, 15..150\}$ to be the expected radii per hough transform. The pseudocode looks like:

   $c_x$ $\leftarrow$ $\emptyset$, $c_y$ $\leftarrow$ $\emptyset$

   for all ($r_i \in \mathbf{r}$) $\{$

      {x, y} = hough_circles($C$, $r_i$)

      $c_x \leftarrow \{c_x \cup x \}$

      $c_y \leftarrow \{c_y \cup y \}$

   $\}$

Depending on the peaks of the accumulator, Hough might report multiple center locations. A simple clustering might help avoiding this (Just average the close-by centers).

The final image with the resulting centers overlayed (in different colors), looks like:

You could customize the radius to be robust to different circles. Moreover, one can use elliptic Hough transform to have a better coverage of the particles (as some of them are not circular). Finally, the much-occluded large circles can be obtained after masking out the small ones - even with the help of a distance transform. If the false circle hypothesis are desired to be eliminated, you can always reduce the Hough sensitivity. Or you could always render the contours of the circle and compare to the detected edges.


aliasing to a 50 Hz sine wave sampled at 5000 Hz sampling rate? I sampled a 50Hz sine wave (its not perfectly a sine wave, it's a data coming from current sensor using Arduino).

Using real time graphing technique by sending the data from Arduino to the PC I managed to reconstruct the signal as attached here

The second experiment I logged my data into an SD card, my sampling rate was 5000 and when I tried to reconstruct the signal in MATLAB, I got this graph which doesn't construct the original signal back, but rather its strange as bellow:

Notice the maximum values for the both graphs.

Why the signal is not correctly reconstructed? I sampled it at 5000 Hz which is almost 10 times the original signal frequency and it meets Nyquist criteria ... And as a result of wrongly constructed signal the FFT is indicating the frequency of zero not 50 Hz.

The MATLAB code

filename = 'C:\Users\Sabir\Desktop\internship Sri Wawasan\testFFT.csv';

% thee file

fid = fopen(filename,'rt');
[data]=textscan(fid, ' %f',...%the values of the sensors as cell array
       'delimiter',',')  ;
fclose(fid);

sensorMatrix=cell2mat(data);

% the values of the sensors as matrix

%taking some point from the sensors matrix

y=sensorMatrix(1:3000,1);% array of one chanel

 plot (y)

L=length(y)

Fs = 5000;                    % Sampling frequency

T = 1/Fs;                     % Sample time              % Length of signal

t = (0:L-1)*T;                % Time vector

plot(Fs*t,y)

and finally here is the FFT, it has 100 Hz and 150 Hz components

You do have zeroth frequency in your signal indeed. If you don't, your signal would be centered around zero volts. Zeroth frequency is a DC component.

If you want to measure the frequency of an AC component, center the signal around zero volts (subtract the mean). By the way, FFT is not a good instrument to measure frequency of an arbitrary signal because it uses a fixed set of harmonics, which may or may not coincide with your signal (and FFT will give you a maximum on a nearest one).


What is the difference between a one dimensional and a two dimensional signal? If a signal depends on only one variable then we call it one dimensional, and if a signal depends on two variable we call it a two dimensional signal.

But when we represent an one dimensional signal, we use two axes(amplitude vs. time).
In case of a two dimensional signal, we use two axes (for example, the x axis and the y axis in an image).

How can I understand the difference between a one and two dimensional signal?

Basically the dimensional number refers to the number of independent variables (input).

In one dimensional signal $f(t)$, the amplitude $y=f(t)$ is the dependent variable (output), and there is only one independent variable $t$.

In two dimensional image $f(x,y)$, the two independent variables are $x$ and $y$, and the dependent variable $f$ is also called intensity. If you want to make it similar with amplitude vs time plot as you do with one dimensional signal, you can plot a intensity vs x vs y scattered dots with three axis as well.


Magnitude of intermodulation distortion products? Assuming I have a set of frequencies $\mathbf F$, and I pick 2 frequencies $f_1$ and $f_2$ randomly from the set, generate sine waves for them, and pass the signals through a non-linear distortion process while adding them.

Is there any way for the harmonics or intermodulation distortion products to have a magnitude larger than the magnitude of any of the two frequencies, when calculating an FFT?
If yes, under what conditions, and how can I avoid that this happens?

Certainly the amplitude of a harmonic can exceed the amplitude of the fundamental,
or the amplitude of an intermodulation product can exceed the amplitudes of the individual signals.   But this occurs when the system is so nonlinear that it would
make no sense to use it as an approximation to a linear system.
Indeed, the very
fact that such a nonlinear system is being used (intentionally)
means that it has some properties
that are useful in the application at hand.

Example: If the system is a squarer that produces output $x^2(t)$, then
with $x(t) = \cos(2\pi f_c t)$, the output is entirely the
second harmonic (and a DC term). Since the fundamental is entirely
absent from the output, the harmonic amplitude is infinitely larger
than the fundamental.  If two sinusoidal inputs of frequencies $f_1$ and
$f_2$ are added together before the
squaring operation, then the output has second harmonics at frequencies
$2f_1$ and $2f_2$, intermodulation distortion products at frequencies
$f_1+f_2$ and $|f_1-f_2|$ and no fundamentals at frequencies $f_1$ and
$f_2$.

Moral: Don't use a squarer and pretend it is a linear addition circuit
or linear amplfier. If you are knowingly using a squarer, it is for
some other useful properties (e.g. frequency doubler) of the circuit that
are important to you right now.


Design FIR filter with half zero coefficions I noticed a lot of hardware IC use special FIR filters with half zero of coefficients to minimize multiplication operations, especially in decimators/interpolators. E.g.
[ 6 0 -19 0 47 0 -100 0 192 0 -342 0 572 0 -914 0 1409 0 -2119 0 3152 0 -4729 0 7420 0 -13334 0 41527 65536 41527 0 -13334 0 7420  0 -4729 0 3152 0 -2119 0 1409 0 -914 0 572 0 -342 0 192 0 -100 0 47 0 -19 0 6]
or
[-12 0 84 0 -336 0 1006 0 -2691 0 10141 16384 10141 0 -2691 0 1006 0 -336 0 84 0 -12]
I took these coefficients from IC datasheet. I see the center coefficient is 0.5 and other odd coefficients are zero. It effectively halves the mul operations.
But I can't find what the type of the FIR filter is. Is there any good paper? How to create a filter with desired number of taps with this feature?

Oh, these are "Nyquist (M) filters", in your case M=2. You'll very often find them in mutlirate systems, especially as L-th-band-filters (i.e. filters that only let through 1/L of the Nyquist bandwidth) in decimators.
This is very fundamental for multirate systems, not really "research paper stuff", so your best bet is hence a good book on Multirate Systems. Fliege has written a classic "Multirate Digital Signal Processing", harris did, too, with "Multirate Signal Processing For Communication Systems" and Vaidyanathan seems to have written a popular book called "Multirate Systems and Filter Banks", but I haven't ever read that.


What are the units of the product of two signals? In a standard mixer, we multiply two signals $x_1(t)$ and $x_2(t)$, each with units of volts $(V)$, and get at the output:

$$y(t)=x_1(t)\,\cdot\,x_2(t)$$

I have always thought that the output of mixing operation was another voltage signal. But how can this be? Looking at the units of the factors, we multiply a quantity of volts, $x_1(t)$, by another quantity of volts, $x_2(t)$. So shouldn't the units of $y(t)$ be $V\cdot V=V^2=W$?

If the multiplier takes two voltages as input and returns a voltage as output, then there is necessarily a constant involved, with units of [1/V].

Take for example, AD633 (which was the first search result). The output is the product of the 2 inputs times a constant:

$V_{out} = \frac{V_1 \times V_2}{10V}$

So the output units are Volts.


about sampling frequency I am trying to make a DFT on a signal with frequency f=50 MHz main component plus some noise.
As far as I know if I sample it at F=100MHz I should be able to get a proper plot of the DFT since F=2f and I should see a peak at 50MHz, this instead does not happen,see first figure.

But instead I obtain what expected when I am leaving unchanged all parameters of DFT (sampling period and number of samples) but halving the signal frequency to f=25MHz, see second figure

Can anyone explain to me why this happens? From Nyquist theorem I should be able to see the 50MHz component if sampling at 100MHz

Thanks for helping

The second plot is correct, 25 MHz tone where you halved your input frequency to 25 MHz. I suspect that you may be sampling synchronously a real signal such that the samples are all landing close to the zero crossings (noticing the dB scale on the horizontal axis on the first plot that there is minimum signal content in your FFT). If this was done with complex sampling the tone would always show up somewhere in the digital domain from 0 to $F_s$. But with real sampling we can get a null at $F_s/2$ as you are seeing with synchronous sampling at or close to the zero crossings.


Recursive form representation of FIR filter Do we have methods to represent FIR filter in recursive form?

One of the methods that I know is frequency sampling method.

Can someone provide me any method?

The frequency sampling representation of a length $N$ FIR filter's transfer function

$$H(z)=\frac{1-z^{-N}}{N}\sum_{k=0}^{N-1}\frac{\tilde{H}[k]}{1-e^{j2\pi k/N}z^{-1}}\tag{1}$$

is a systematic way of finding a recursive implementation for any FIR filter. The coefficients $\tilde{H}[k]$ in $(1)$ are the DFT coefficients of the impulse response $h[n]$. I do not know of any other such systematic way.

However, there are special cases for which a recursive implementation can be found by using standard finite sum identities. One such class of finite impulse response transfer functions is

$$H(z)=K\sum_{n=0}^{N-1}a^nz^{-n}=K\frac{1-a^Nz^{-N}}{1-az^{-1}}\tag{2}$$

A special case of $(2)$ is the moving average filter

$$H(z)=\frac{1}{N}\sum_{n=0}^{N-1}z^{-n}=\frac{1}{N}\frac{1-z^{-N}}{1-z^{-1}}\tag{3}$$