Can all pass filers have a pole/zero at the origin of the root locus I came across this root locus diagram in my DSP class. It has a pole at the origin, another pole within the unit circle and a zero just outside the unit circle. I understand that the poles and zeros should be complex conjugates for an all pass filter. So will a pole or a zero at the origin prevent the system from being an all pass?My professor tells me that the pole at the origin is allowed. So that is basically leaving me confused.Can anyone explain how this is possible?

a delay is also an APF of sorts.  it passes the same magnitude and it affects the phase (in a linear manner w.r.t. frequency).  a delay of one sample has transfer function

$$ H(z) = z^{-1} = \frac{1}{z} = \frac{1}{z-0} $$

which means there is a pole at $z=0$.  and it is true that $\left|H\left( e^{j \omega} \right)\right| = 1$, so it's an APF.

so, if the remaining pole/zero pairs are at reciprocal values this extra delay changes the phase of your APF a little, but not the magnitude response, so i think, strictly speaking, that it's still an APF.


How to estimate filter coefficients of an all-pole system? Suppose we have measured the frequency response of a system that is known to be all-pole; measuring impulse response is not possible. What are the methods available, if any, to estimate the coefficients of the underlying IIR filter?

EDIT: Frequency response exists only as measured data, so closed mathematical expression is not available, nor is the phase information.

Half of the information about a system is in the phase response, and half of it in the amplitude response. There's no way to reconstruct a filter if you know only either half, unless you strongly restrict the poles mathematically.


Calculation of an impulse response of h[n] I am currently looking at the z-transform and am using a great youtube reference to help me, however I am struggling on some basic step. How do I get the impulse response array of h[n] = [ ... ] shown in the lecture at the following URL - https://www.youtube.com/watch?v=dEJp46SFgV4&ab_channel=DavidDorran

I can do the following which I believe to be correct but do not know how to get the array out as shown in purple.

That purple array is giving the impulse response over time. You can get it directly from the difference equation. Assume initial rest, $y[-1]=0$, then write out the impulse response for $n=0, 1, 2, ...$. If you do that you will get:
\begin{align}
 h[0]&=1+0 \\
 h[1] &= 0 + 0.5 \\
 h[2] &= 0 + 0.25 \\
 h[3] &= 0 + 0.125 \\
 ...
\end{align}


Estimation of Input Signal to Obtain the Desired Output Signal for an Unknown Filter Suppose $h(n)$ is a finite impulse response which is unknown. We can feed any input signal $x(n)$ into the system and observe the corresponding output signal $y(n)$.

From this, is it possible to estimate the input signal $x^*(n)$ which when passed through the filter gives the desired output signal $y^*(n)$ ?

I think a Least Mean Squares filter could be used here. Could anyone point out how one proceeds doing this ?

Thanks !

In the most simple case, just to give intuition about the problem, it is really easy.

In the Frequency Domain:

$$ {Y}^{\ast} \left( \omega \right) = H \left( \omega \right) {X}^{\ast} \left( \omega \right) \Rightarrow {X}^{\ast} \left( \omega \right) = \frac{ {Y}^{\ast} \left( \omega \right) }{ H \left( \omega \right) } $$

Since $ {Y}^{\ast} \left( \omega \right) $ is known all needed is $ H \left( \omega \right) $.

Yet since we have access to a black box of $ H \left( \omega \right) $ we can set input of a known signal $ X \left( \omega \right) $ and have the output $ Y \left( \omega \right) $ which will give us $ H \left( \omega \right) $:

$$ Y \left( \omega \right) = H \left( \omega \right) X \left( \omega \right) \Rightarrow H \left( \omega \right) = \frac{ Y \left( \omega \right) }{ X \left( \omega \right) } $$

As others pointed out it is called Deconvolution and as you mentioned it can be done using Least Mean Square (LMS) Filter.


MATLAB IIRs -- extracting and scaling the coefficients of second order sections suppose I generate the following filter:

Fs = 32000;  % Sampling Frequency
N     = 4;    % Order
Fpass = 400;  % Passband Frequency
Apass = 1;    % Passband Ripple (dB)
h  = fdesign.lowpass('N,Fp,Ap', N, Fpass, Apass, Fs);
Hd = design(h, 'cheby1', 'FilterStructure', 'df2sos');

"Hd" at the interactive prompt gives me:

 sosMatrix: [1 2 1 1 -1.97232454484577 0.978340933957857;1 2 1 1 -1.94676435093801 0.948444394718458]
 ScaleValues: [0.00150409727802169;0.000420010945111697;0.891250938133746]

I notice that freqz(Hd.sosMatrix) gives me a wildly off magnitude estimation (but with the correct frequency cutoff).

My question is, how do I use this Hd.ScaleValues factor? There are 3 values for two sections.

These scale values are applied before and after each section, so for $N$ sections you have $N+1$ scale values. Check the Matlab documentation, there must be some related property that you can set to 'scalar', and this will give you a single total scaling value. In any case, the total scaling is just the product of the three values in your case. So you have to multiply the b coefficients of one filter section by the product of all scale values, then the total scaling will be OK, i.e. the maximum frequency response magnitude will be 0 dB.


How matching pusuit algorithm works I could infer from books and wiki that a matching pursuit finds signal approximations from a redundant dictionary. Let f be the signal, ai be the basis and Ri be the residue.. here we choose the inner product of f and ai to be max such that Ri is less. What do we do next. what do we iterate?

The point is that the dictionary is redundant or overcomplete, i.e. there is no unique expansion of a given signal, unlike the expansion in an orthonormal basis. This means that one tries to find "the best" expansion of a given signal, e.g. such that the error between the signal and its expansion is small while most expansion coefficients are close to zero, which is advantageous for signal compression. The matching pursuit algorithm is an iterative algorithm trying to find an approximation of the best expansion of a signal given a dictionary. Here you can find another useful overview of the matching pursuit algorithm.


Detecting a curve with certain characteristics in a spectrogram I have a spectrogram built by breaking a signal into chunks and taking the FFT of these chunks. An example spectrogram is:

In this spectrogram I am looking for a curve (aka pitch contour) that:

Starts at zero
Increases to an unknown value
Stays at that unknown value for at least one minute
Decreases back to zero

It's fairly easy to program this curve detection in the spectrogram above since the curve of interest is brighter than other potential curves.  However, this brightness difference won't exist in a noisier data set.

How can I use the constraints I know about the curve (starts at zero, increases, ...) to reliably find the curve without relying on the brightness?

This answer might sound hard initially, but once understood it is intuitive and certainly deals well with the noise.

One smart way to model such constraints is either through HMM's or through MRFs. The design is similar, but for HMM I would give you some hint:

Prepare your transition matrix in such a way that Viterbi decoding would provide you the resulting curve. For example, the first state will be at 0 with certain variance, since it starts there. Then it is only allowed to increase, within a given tolerance. This tolerance would be obtained from a Gaussian prior. Say, the slope of the line (thus the variance of the Gaussian) would range from 1 to 5. Finally, you should only allow for a decrease or keeping constant. Similarly, for every other constraint you could do this. Here allowance means that the probability of transition from the constant state to an increasing state is 0. Otherwise you have a probability (might be due to any distribution, prior etc.)

You could also make use of Hierarchical HMMs to elaborate your model with such tracking parameters:

http://en.wikipedia.org/wiki/Hierarchical_hidden_Markov_model
http://www.cs.princeton.edu/courses/archive/spr06/cos598C/papers/FineSingerTishby1998.pdf


Implementing FFT in Matlab I am facing a trouble implementing the fft and appreciate your advice. My data structured as follow:
The data are retrieved from the share market for one day change in the share price. The time range from (11:00:00 to 15:00:00) this means 16200 Sec. In matlab its represented as serial so it will be from 0.4584 to 0.6458. Data are not sampled in fixed freqency. Price is sampled only when there is a change in price. For example ( it can be sampled in each second and sometimes in 30 Seconds interval).

For calculating the frequency spectrum of an non-uniformly sampled signal I see two options:

Interpolate the data in order to obtain evenly spaced samples before taking the FFT as suggested by Hasan. Use the Matlab function interp1. If you have the Data Analysis toolbox, the resample function of the timeseries class might be worth a look.
Use the Lomb–Scargle periodogram. There's no builtin Matlab function but an implementation from Matlab File Exchange exists.


Background averaging not working I'm working on a problem where I'm trying to use background averaging to see what the frame looks like in a beehive behind the bees. I've been trying to use standard averaging (I take an image every 1 or 4 seconds over the course of my 3 hour video) as well as other techniques like mog and mog2. The problem is that there is an extremely high density of bees in the centre that is tightly packed and barely moving, so while I end up with a great view of the background around the middle, the centre turns into a bit of a smear like this:

Does anyone have any advice on what I could do to improve things? I spoke to a researcher who suggested something called "sparse coding" as one possible options and I'm curious if OpenCV has anything like that or if anyone knows of any other techniques I could try.

I would try techniques related to low-rank approximations of your data $X$. For instance, the paper GoDec: Randomized Low-rank & Sparse Matrix Decomposition in Noisy Case attempts to decompose
$$X = L + S + G$$
with a low-rank part ($\mathrm{rank}(L) \le r$), a sparse part ($\mathrm{card}(S) \le k$) and a residual $G$.

It has been employed for background modeling, with backgrounds related to the low-rank part, plus sparse en independent objects.

Some video demos are:

GoDec and Robust PCA for background modeling on Escalator
GoDec and Robust PCA for background modeling on Hall

There are alternatives in the field of robust matrix factorization, and I understand it as an instance of sparse coding.


Effect of sign change on ROC & z transform? When we take ztransform of unit step it is z/[z-1]
And our ROC is |z|>1
But if some how the minus sign between z and 1 changes to +, will our ROC still be same as old ( |z| >1 ) or will it be opposite?

Secondly the ztransform and ROC of u(-n) and u(n) are same or not?

You need to understand that a function such as $F(z)=\frac{z}{z+1}$ doesn't have a region of convergence (ROC). Only an infinite series has a ROC. E.g., the $\mathcal{Z}$-transform of a sequence $x[n]$ can be expressed by the series

$$\sum_{n=-\infty}^{\infty}x[n]z^{-n}\tag{1}$$

and this expression only makes sense if the series converges for certain values of $z$. The ROC is the region of the $z$-plane for which the series $(1)$ converges.

If you turn the problem around and if you're given a function of $z$ like $F(z)$ defined above, you can ask yourself if $F(z)$ is a valid $\mathcal{Z}$-transform of some sequence $f[n]$. It turns out that $F(z)$ is the $\mathcal{Z}$-transform of two different sequences, but with different ROCs. So a $\mathcal{Z}$-transform is only completely defined by the expression for the function $and$ the corresponding ROC.

Coming back to the function $F(z)$, it's straightforward to show that it is the $\mathcal{Z}$-transform of the sequence $f_1[n]=(-1)^nu[n]$ with the ROC $|z|>1$, and of the sequence $f_2[n]=-(-1)^nu[-n-1]$ with the ROC $|z|<1$.

Concerning your last question, as mentioned above, the expression for the $\mathcal{Z}$-transform including the ROC is unique, so the two sequences $u[n]$ and $u[-n]$ cannot have the same $\mathcal{Z}$-transform and ROC. It's actually pretty easy to come up with the $\mathcal{Z}$-transform of $u[-n]$:

$$\begin{align}\mathcal{Z}\big\{u[-n]\big\}&=\sum_{n=-\infty}^0z^{-n}\\&=\sum_{n=0}^{\infty}z^n\\&=\frac{1}{1-z},\qquad |z|<1\end{align}$$

which is of course different from the $\mathcal{Z}$-transform of $u[n]$. Note that the algebraic expression as well as the ROC are different.


Is there any good book regarding 3D reconstruction from 2D images? I prefer to work in MATLAB I am working on a project regarding 3D reconstruction of 2D images. I have understanding of finding and Mathcing the keypoints and its algorithms; but i am not able to replicate the theory in MATLAB code. Is there any Book or paper which is basic for this problem?

There are a lot of good books. Multiple View Geometry in Computer Vision,

"Computer Vision: Algorithms and Applications", but anyway, replication of the theory is the necessary stage to me.

You can start from the finishing course "3D Computer Vision" from the CTU. They published lectures and labs - in MatLAB.

P.S. Here is a ready implementation http://www.fit3d.info/.


Doppler Model, how does it affect the fading impulse response? I have been reading a paper that mentions Doppler model that may be included in the fading characteristics of impulse response. One example of bell shaped Doppler Spectrum is the following

$$S(f) = \frac{1}{1+A\frac{f}{f_d}}\,\,\,\,\,\,\,\, |f|\leq f_{\max}$$

where $f_d$ is the Doppler Spread, $A$ is a constant, $f_{max}$ is maximum frequency component of the Doppler power spectrum. Does anyone know how the following function affects the channel model?

Thanks

This is not a full answer, but I'll try to provide some intuition. First, imagine a static channel -- neither the antennas nor the reflectors around them move. In such a channel, the impulse response is always the same -- it's time invariant. The channel characteristics are fixed.

Now, introduce movement of the antennas and/or the reflectors. The channel impulse response now varies with time. The channel coherence time indicates how fast the channel impulse response changes, and for how long can we consider the channel to remain relatively static.

If there is movement of the antennas and/or the reflectors, we also get a Doppler shift in the transmitted signal. The more movement there is, the higher the shift.

The channel coherence time and the Doppler shift are, then, two manifestations of the same phenomenon -- movement. If one increases, the other decreases. The easiest way to introduce this into the channel model is to assume that the channel is fixed for a certain amount of time, after which it changes to a completely different, independent channel.


What does it mean for a signal to alias to another? I know the definition of this is simple and I do understand the concept of aliasing (I think) in the way that multiple signals can be aliases of each other. For example, I get it when some one says: sampling of a 7-kHz sinewave at a sample rate of 6 kHz will provide a discrete sequence of numbers whose spectrum ambiguously represents tones at 1 kHz, 7 kHz, 13 kHz, 19 kHz, etc. However, how do you determine what aliases to what? Obviously I have a terminology problem. I understand that signals can be aliases of each other (right??), but cannot see what signal aliases to what other signal.

Sampling of a signal leads to repetition in the frequency domain. Typically, the frequency domain is considered from $-\pi/2$ up to $\pi/2$ or from $0$ to $\pi$, where $\pi$ corresponds to the sampling frequency.

For a real-valued signal, the spectrum is also symmetric around $n\pi/2$ for any integer $n$. For your example where a sine wave with frequency $7$kHz is sampled at $6$ kHz, this means that the sampled signal contains the frequencies $1+6\cdot n$ kHz and the symmetric/mirrored frequencies $5+6\cdot n$ kHz, e.g., 1 kHz, 7 kHz, 13 kHz and 5 kHz, 11 kHz (and also the negative frequency parts.

You speak of 'aliasing over another signal' when the aliases overlap with other signal content. For example, when you have a signal consisting of the sine wave of 2 kHz and 7 kHz, which you sample with a rate of 6 kHz, then you observe the frequencies 2 kHz; 4 kHz and 1 kHz; 5 kHz. These frequencies do not overlap. When the original signal contains 1 kHz and 7 kHz and you sample at a rate of 6 kHz, then the '7 kHz signal aliases over the 1 kHz signal'.

Of course this principle works also with other signals, such as narrowband, low passed and high passed signals.


1-D discrete cosine transform(DCT) This was a past year paper question. Not sure how to answer it.

Question: The 1-D discrete cosine transform(DCT) of a sequence f(x), x =0,1,...,N-1 is

F(u) = c(u)* Summation of f(x) cos ( (2x+1)u*pi/2N)

show that 1-D DCT of the sequence

g(x) = f(N-1-x), x = 0,1,...,N-1

can be expressed as G(u) = (-1)*F(u), u=0,1..., N-1

My solution (If wrong please correct me)

G(u) = c(u)* Summation of f(N-1-x) cos ((2(N-1-x)+1)u*pi/2N)
     = c(u)* Summation of f(N-1-x) cos (2N-2x-1)u*pi/2N)

If $g(x)= f(N-1-x)$ then
$$\begin{aligned}
  (\operatorname{DCT}\:g)(u)
    =& c(u)\cdot\sum_{x=0}^{N-1}g(x)\cdot\cos\bigl((2x+1)u\cdot\tfrac\pi{2N}\bigr)
  \\=& c(u)\cdot\sum_{x=0}^{N-1}f(N-1-x)
        \cdot\cos\bigl((2N-2N+1+2x)u\cdot\tfrac\pi{2N}\bigr)
  \\=& c(u)\cdot\sum_{(N-1-x)=(N-1)}^{0}f(N-1-x)
        \cdot\cos\Bigl(\bigl(2N-1-2(N-1-x)\bigr)u\cdot\tfrac\pi{2N}\Bigr)
\end{aligned}$$
substitute now $N-1-x=:\xi$,
$$\begin{aligned}
  (\operatorname{DCT}\:g)(u)
    =& c(u)\cdot\sum_{\xi=0}^{N-1}f(\xi)
        \cdot\cos\Bigl(\bigl(2N-1-2\xi\bigr)u\cdot\tfrac\pi{2N}\Bigr)
  \\=& c(u)\cdot\sum_{\xi=0}^{N-1}f(\xi)
        \cdot\cos\Bigl(\pi\bigl(\tfrac{2N}{2N}-\tfrac{2\xi+1}{2N}\bigr)u\Bigr).
\end{aligned}$$
In the cosine, we have now a shift by an integer multiple ($u$) of $\pi$. Such a shift in a sine or cosine function is equivalent to multiplying the function with $(-1)^u$, So the result is
$$\begin{aligned}
  (\operatorname{DCT}\:g)(u)
    =& c(u)\cdot\sum_{\xi=0}^{N-1}f(\xi)
        \cdot\cos\bigl(\pi\tfrac{2\xi+1}{2N}u\bigr)\cdot(-1)^u
  \\=& (-1)^u\cdot(\operatorname{DCT}\:f)(u).
\end{aligned}$$
Not quite the expression you wanted to show, but this should be the correct one.


Show that $S[n]= \frac1M \sum e^{\frac{2\pi l n}M}$ for $l = 0$ to $M - 1$ I am learning about image processing. Currently, I am reading about upsampling and down sampling.
Here is some background what I am reading:

I was wondering if someone show me a proof of how equation (70) is true. It simply says "It can be shown that".

Thanks!

HINT:

First of all, Eq. $(70)$ is wrong, there's an imaginary unit $j$ missing in the exponent. If you add that missing $j$, note that for $n=kM$, $k\in\mathbb Z$, you have $e^{j2\pi ln/M}=1$, so the result of the scaled sum equals $1$. In order to see that the sum equals zero for all other values of $n$, use the formula of the geometric series.

This answer is also relevant to your question.


Confusion in proof of Inverse Laplace Transform For the proof of inverse Laplace transform, we change the integral from $\omega$ to $s$. I want to know the reason why we need to change the integral?

To summarize the discussion:

The usual substitution is $s = \sigma + j \omega$ where $\sigma$ is the real part of the $s$ variable and $\omega$ is the imaginary part.
The equation in the image is for the Fourier transform, not the Laplace transform. The Fourier transform can be thought of as the Laplace transform evaluated on the imaginary axis ($\sigma = 0$).
The differential $ds$, when looking at real and imaginary parts distinctly, becomes $d\sigma + j d\omega$.
Any differential is an infinitesimal (very small) change in that variable. $dx$ is a small change in $x$.


Optimal $ n $ -th Order IIR /AR Approximation of a Moving Average Filter I would like to approximate a moving average filter with an IIR filter of much lower order than the tap-length of the moving average filter. Optimality shall refer to the $L_2$ norm of the impulse response (up to changes in the total gain). Ideally, I'd also like to constrain the impulse response to be monotonic (i.e. no ripple), but I fear that this requirement along will already fully constrain the filter. If that is the case, I might be willing to take some ripple in turn for a sharper roll-off.

(Update: I just realized that the first part of this covering CIC structures is basically what Hilmar has already answered-- I'll leave this up since it offers more graphics and details in case that his helpful to anyone but it is indeed the same answer)

This may not be optimum and although a feedback structure is involved it is strictly not IIR, but want to point out the efficiency of the cascade integrator-comb structures for implementing boxcar moving average filters- and importantly how the results match exactly and not an approximation. So a 600 tap moving average filter can be implemented with simply one subtraction after a 601 sample delay followed by an accumulator.

These structures are quite popular when resampling is involved (often after a moving average the higher sampling rate of the input is no longer required and a lower output sampling rate would be even more efficient).

The equivalency is easily proven from the Z-transform for each of these structures. Specific to the example structures in the graphic:

The Z-transform for the accumulator is:

$$\frac{1}{1-z^{-1}}$$

The Z-transform for the comb (delay and subtract) filter is:

$$1-z^{-4}$$

And the Z-Transform for the moving average filter is:

$$1+z^{-1}+z^{-2}+z^{-3}$$

Combining confirms the equivalence:

$$1+z^{-1}+z^{-2}+z^{-3} = \frac{1}{1-z^{-1}}(1-z^{-4})$$
$$(1-z^{-1})(1+z^{-1}+z^{-2}+z^{-3})= (1-z^{-4})$$

Another favorite of mine as a true IIR alternative to a moving average filter is the exponentially weighted moving average filter, where assuming sufficient bit precision, very tight (in other words effectively long delay) filtering can be achieved by using an alpha that approaches 1. This is demonstrated in the following graphics:

In fixed precision structures the order of the two operations shown in the block diagram of the exponentially weighted moving average filter is very important! In most cases when dealing with fixed precision it is best to allow digital filters to grow the signal first and do the scaling down after the accumulation. This also explains the motivation for extended precision accumulators- the accumulator shown needs to be designed with sufficient precision to prevent any overflow for the given $\alpha$ value chosen (the Infinite Summation identity $\sum_{n=0}^\infty \alpha^n =\frac{1}{1-\alpha}$ is useful in this case to determine worst case maximum value; resulting in $gain = \frac{1}{1-0.99}$ or 100 in the example shown, so in this case we would need to add $log2(100) = 7$ bits in the accumulator). Also clear from the summation identity is how the scaling term itself was determined such that the gain of the structure is 1.
The order of these operations is no different than the reason to put an amplifier prior to an attenuator in an analog system to get the best noise figure, while paying attention to clipping and maximum input signal, for those familiar with that.


continuous chain of ones or zeros for transmit If my data has long chain of 1's or 0's. Can the data be modulated and transmitted as it is. Will this have any effect of  the receiver adc as this may increase dc value of the signal or any other block in the receiver chain?
should the data be scrambled?

A dominant concern with long runs of ones or zeros is the ability of the receiver to establish symbol timing. The receiver clock is independent from the transmitter clock, so without transitions in the data the receiver will not be able to determine the clock used in the transmitter for each symbol, and will ultimately drift off. Whitening the data has other benefits for equalization and overall efficiency.


How to detect a long signal inside an even longer one? I'm trying to detect a long pattern inside a long recorded audio. In order to do so I am implementing a matched filter in time by performing a FIR filtering using the coefficients of my pattern inverted in time as filter coefficients and doing all the processing in frequency domain. I'm ensuring periodic convolution is equivalent to linear convolution by carefully selecting the FFT size named NFFT so that there's no time aliasing at the output.

$ ifft(fft(patternCoeffTimeInverted,NFFT).*fft(Segment,NFFT),NFFT)$

I am familiar with the existence of overlap-add and overlap-save algorithms to optimally compute the FFTs I would use to implement the convolution between pattern and audio. The problem I'm facing is that normally these algorithms operate with a manageable number of filtering coefficients A and the large audio is divided into chunks of size B for easier processing. My question is what can I do if both the filtering coefficients and the audio are large? Would I have to nest overlap-adds by taking the first chunk of audio B and use it as coefficients A' in a second overlap-add while dividing the original coefficients from A into chunks just to have an output to be used in an iteration of the first overlap-add and then do this for every iteration of the first overlap-add?

Using the matched filter is right approach and typically a full length filter should be used.

One of the advantages is that  filter gain is approximately proportional to $10\log(T)$.

If you don’t need full gain you can use a chunk of your matched filter and still have acceptable detection performance.

One technique in active sonar is to use a broken up matched filters because returns tend to be broken up. A target isn’t a point reflector.

You can also use linearity and time invariance to break up the long convolution. Say take the first half of your matched filter and convolve it with your long time series. Then do the same with the second half matched filter. If you align the 2 convolved time series, offset by the length of the first half filter, adding the convolved time series is equivalent to convolution by the full matched filter


Implementation of the filter function (Matlab) C/C++ I am trying to apply some filtering in C/C++ code, I have designed an algorithm using Matlab and it works well, so now I am moving it to native code.

I am interested in the filter function.
In my matlab script it is used as follows.

[y_blk, fir1_state] = filter( h, 1, x.*c, fir1_state );

Where the h argument is coefficients generated by firpm function.

I have tried to find any implementation, but unfortunately no results.

I would be very grateful if someone could provide an implementation of this function

Or maybe there is an existing implementation of this function in some library.

Thanks

FIR filtering is an operation whose implementation hugely benefits from utilization of architecture-specific instructions, compiler intrinsics, etc. Hence, if your only goal is to get a performant FIR algorithm, you should try to find a dsp library specifically for your target architecture. Those are most likely written in assembly language with a C interface and you can expect an orders-of-magnitude better performance than you would achieve with a generic C implementation of a transversal FIR.

If the latter is exactly what you're looking for - a generic FIR implementation in C - let me point you towards https://iowegian.com/. They offer a commercial filter design software called ScopeFIR. You can download and install a free trial version of it (no registration). In the program directory, you can find a source file "FirAlgs.c" that offers a variety of well-documented FIR implementations. It's written for MSVC but should be portable to basically anywhere without hassle.


Pinhole camera model from houdini parameters From a simulation with Houdini software I retrieved these parameters:

Camera position: -0.675839, 33.5945, 0.0318854
Camera Rotation: 0.0318854, 92.4693, 0.1
Focal Length: 52.5172 mm
Image size: 1920 x 1080
Pixel Aspect Ratio: 1
x aperture = 41.4214 mm

I am not sure about how to build the system for transforming 3d coordinates of the scene to 2d image coordinates, in particular how to derive and combine rotation-translation matrix, camera-matrix transformation, and image transformation.

Do you have any suggestion?
Thank you

This is a long topic to fully explain. I will try to write shortly, so please excuse the brevity.

Standard computer vision projection (ignoring distortion like Houdini) follows:

$$
\mathbf{x} = \lambda \mathbf{K}[\mathbf{R}\mathbf{X} +\mathbf{t} ]
$$

$\mathbf{R}$ is a $3x3$ orthogonal matrix, $\mathbf{t}$ is a $3x1$ translation vector. Camera position $\mathbf{C}$ is given by $\mathbf{C}=-\mathbf{R}\mathbf{t}$. $\lambda = x_z$ is a scale factor used for de-homogenization. The projection can also be compactified into a $4x3$ matrix:
$$
\begin{align}
\mathbf{x} &= \lambda\mathbf{P}[\mathbf{X}^T 1]^T\\
\mathbf{P} &= \mathbf{K}[\mathbf{R} | \mathbf{t}]
\end{align}
$$

Rotation is sometimes represented as a Rodrigues vector. To simply convert from one representation to another you might use cv::Rodrigues. Sometimes, (for some unknown reason!) people can also choose Euler angles for rotation representation. You should check which one Houdini uses, hopefully Rodrigues - if not you should compose $\mathbf{R}$ from 3 angles. To obtain $\mathbf{t}$ explicitly, just plug in $\mathbf{C}$ and $\mathbf{R}$ into $\mathbf{C}=-\mathbf{R}\mathbf{t}$ to solve for $\mathbf{t}$.

The only remaining unknown is then the $3x3$ camera instrinsic matrix, $\mathbf{K}$. In computer vision:

$$
\mathbf{K} =
\begin{bmatrix}
f_x & s & c_x \\
0 & f_y & c_y \\
0 & 0 & 1 \\
\end{bmatrix}
$$

The principal point $(c_x, c_y)$ is sometimes initialized to be $(\frac{w}{2},\frac{h}{2})$, which you can obtain from image size. Usually $s=0$. $f$ is just a bit more complicated:

Standard vision libraries (e.g. OpenCV) store the focal length in pixels, as in $\mathbf{K}$. In general, the relation to physical focal length is:
$$
f_{in-pixels} = \frac{w * f_{mm}}{w_{CCD-in-mm}}
$$
where $w_{CCD-in-mm}$ is the physical sensor size. In your case, I suspect that the aperture is $w_{CCD-in-mm}$, and focal length is the same as $f_{mm}$. If you plug these values you will get the focal length in pixels: $f=f_{in-pixels}$. Not sure about your units though.

Moreover, I assume that with $1$ x $aperture$ $= 41.4214$ you refer to a camera with square pixels with aspect ratio 1. Therefore it will not affect if you do the same calculation using height instead of width: $f_x=f_y=f$.

I hope this clarifies a bit. Let's keep refining this answer if you provide more info.


Bandpass with large group delay I read that a simple delay filter that delays all frequencies requires at least $\frac {N-1} {2f_s}$ taps. In my case this would result in way too many taps (20000).
Can I achieve this with significantly fewer taps (maybe around 1000, smaller by factor 20) if I either:

only want to distortionless-free delay a bandwidth of $f_s/100$ and don't care about the filter response outside of that band

or is there a way to have the center 18000 taps zero (only 1000 taps on the left and right have coefficients)?

Do you have any tips or could point me to the correction literature for this very specific case? Or are there even python functions that allow me to specifiy these requirements?

This depends a lot on whether your constraints on the number of taps is computation or memory.
If you have a enough memory, you can simply use an FIR with all taps zero except for the last one (Which is $1$). That's equivalent to using a ring buffer and no filter at all.
If you don't have enough memeory, things are more tricky. YOu can get delays that are constant with frequency with parallel/cascaded alpass filters and you can use "schroeder allpass" filters to crank up the group delay effificently. The best way to do this depends a lot on your specific constraints: how much memory do you have, how wide is your frequency range and how "constant" does that delay needs to be.


Sub-Band Extraction Techniques for Complex (I/Q) Signal There is an application which processes Complex I/Q signal of 100 MHz bandwidth and it is required to extract independent sub-bands from this bandwidth at different center frequencies and bandwidths. I wanted to ask what are the techniques which can be used to extract this sub-band data?
Edit: FFT can also be considered as a filter bank. Can the FFT information be used to create independent sub-band data using IFFT etc.?

I wanted to ask what are the techniques which can be used to extract this sub-band data?

That's a description of a filter bank.
Yes, the FFT can be used for such applications. You'll find that OFDM, which powers DVB-T, 4G/5G, WiFi, … (basically all high-speed wireless terrestrial links) does exactly that.
You'll also find that if you find the inherent sinc-shape of subcarriers in an FFT filterbank insufficient, but still need regularly spaced channels, polyphase filterbanks are the way to go.


Validating cross corrolation between sine and cos, shouldn't pahse lag be pi/2? I am trying to use cross correlation instead of FFT to find lags between 2 signals. I am trying to understand cross correlation by using it on sin and cos and expecting phase lag to be pi/2 but it is not:

clc
clear

frequency =1000;
t= linspace(-1000,0.1,1000);

sine_wave = sin(2*pi*frequency*t);
cos_wave  = cos(2*pi*frequency*t);

[xc,lags] = xcorr(cos_wave,sine_wave);
[~,I] = max(abs(xc));

figure
stem(lags,xc,'filled')
legend(sprintf('Maximum at lag %d',lags(I)))
title('Sample Cross-Correlation Sequence')

I am getting a maximum at lag 22 (if you copy and paste the code provided you ll see it in the plot) and I understand that phase lag = 2*pif (time delay) so in my case:

phase lag = 2*pi*1000*22 which is by no means close to pi/2...
What am I missing here?
Thank you

To start with there is a local maxima at 2 lags. If you zoom the x-axis (say axis[-200 200 -500 500]) you can see it more clearly. However, it is not the global maxima and the one at lag 22 is higher.

This is due to the digitization of your sample. If you plot one of your waves against t and zoom the x-axis to say a 200 width box you will see that the peaks clearly have different heights. You therefore get a ringing effect in the resulting cross-correlation.

You could reduce this effect by significantly increasing your sampling rate or decreasing the signal frequency. You currently have ~5 points per $\pi$ of phase which is very low.

You could also adjust your sampling rate/frequency such that an integer number of samples fit exactly in $2\pi$ of phase.


Analytical expression for the eigenvectors of a 3x3 real, symmetric matrix? I am writing an algorithm that process 3D images based on the local moment of inertia.

I have a 3x3 real symmetric matrix, from which I need to find the eigenvalues. I have found a variety of generic algorithm for the diagonalization of matrices out there, but I could not get to know if there exists an analytical expression for the 3 eigenvctors of such a matrix.

Would someone proficient in maths know that?

EDIT

For the record here is what I have found on the question myself.
As Matthias Odisio said, you can't get down to a simple analytical expression as soon as you have a 3x3 matrix.

I have found however a dedicated paper for the special case a 3x3 hermitian matrices, where various numerical specialized approaches are compared:

http://arxiv.org/abs/physics/0610206

Here is the C and Fortran code of the paper:

http://www.mpi-hd.mpg.de/personalhomes/globes/3x3/index.html

Wolfram|Alpha has made such things easy:


ROC curve with classification problems I am working with face anti-spoofing. For classification, the Euclidean distance measure is used. The test image vector is compared with all real image vectors , and the minimum distance is chosen (min.1) . The same vector is compared with all fake image vectors and the minimum value is chosen (min. 2). The minimum rule is applied again to choose the minimum of both min.1 and min.2 to identify the correct class of the test image. My problem is in plotting the ROC curve. I do not have a threshold which can be changed. I am using hard decision rules. How can I plot the curve in this case

Yes. You can.
ROC curve is a Pd vs Pfa plot where Pd and Pfa are the probability of detection and probability of false alarm respectively.


How to check if h(n) is causal, stable? I have $x(n)$ = {$4,-1,-3,1$} and $h(n)$ = {$2,1,3,5$}.

I would like to know how can I check whether

$h(n)$ is stable filter or not and
$h(n)$ is causal filter or not?

Firstly, you do not need to specify an input $x(n)$ to know whether the system is stable or causal, these are properties of the system itslef, I will interpret $h(n)$ you mention as a array describing the system response.

For causality in time domain perspectivefor an LTI system check if $h(n) =0 $ for all $n < 0$, that is the $h(n)$ is essentially a right sided sequence.

For causality in Z domain perspective equivalently, the Region of Convergence is the right side of the outermost pole (outermost in terms of the highest magnitude pole in z domain)

For stability in time domain perspective, the system response $h(n)$ should be absolutely summable, since you have only finite coefficients, this system is stable, any FIR filter is always stable.

For stability from Z domain perspective, the region of Convergence must contain the unit circle.

For stable and causal system all these conditions hold at the same time.


what is an alpha filter? Currently, I'm working on adaptive beamforming using LMS approach, so they change the value of the step factor adaptively in which one of the steps is to pass the weight vector through an alpha filter. SO I'm unable to get any appropriate texts on what it is. Can anyone help me out? Thanks in advance.

My guess is that it is an alpha filter, as defined in the context of alpha/beta filtering.


MATLAB: How to copy column vectors from a matrix into a cell array? I have some sensor data stored in a 1000 x 5 matrix. I'd like to copy each column into a 5 x 1 cell array such that each cell contains a 1000 x 1 vector. I tried the num2cell(mymatrixname,1) function, but this created a cell array of 1 x 5, not 5 x 1. Using MATLAB, how can I do this in a compact way - without using a for loop?
Thanks,
kilosaw


Example of non-equivalence of the two PSD definitions According to the book Introduction to Spectral Analysis by P. Stoica and R. Moses, the power spectral density (PSD) $P(\omega)$ can either be defined as the discrete-time Fourier transform (DTFT) of the covariance sequence $r(k)$, i.e.,
\begin{align}
P(\omega)=\sum_{k=-\infty}^{\infty}r(k)e^{-j\omega k},
\end{align}
or alternatively as
\begin{align}
    P(\omega)=\lim_{N\to \infty}E\Bigg\{\frac{1}{N}\Bigg|\sum_{n=0}^{N-1}x(n)e^{-j\omega n}\Bigg|^2\Bigg\},
\end{align}
and these two definitions are equivalent under the assumption that
\begin{align}
    \lim_{N\to \infty}\frac{1}{N}\sum_{k=-(N-1)}^{N-1}|k||r(k)|=0.
\end{align}
What is an example of a signal for which these two definitions yield different results, i.e., for which the above assumption does not hold?

Work in progress:  wait till I am done before reading (and throwing brickbats!)

This question is difficult to answer without getting into a lot of details about basic signal analysis and Fourier transform theory.

Because of the way my brain works, I will discuss only real-valued continuous-time deterministic signals, and will get into the stochastic and discrete-time stuff later.  The classical Fourier transform theory considers a finite-energy signal $x(t)$ and defines its Fourier transform as
$$X(f) = \int_{-\infty}^\infty x(t) \exp(-j2\pi ft)\, \mathrm dt.\tag{1}$$
Note that $x(t)$ necessarily decays away to $0$ as $t \to \pm\infty$; absent this property, the signal cannot have finite energy as we have assumed.

The energy spectral density of $x(t)$ is defined to be $|X(f)|^2$ which happens to be the Fourier transform of the autocorrelation function $$r_x(t) = \int_{-\infty}^\infty x(\tau)x(\tau+t)\, \mathrm d\tau.\tag{2}$$ Note that $r_x(0)$ equals the finite energy of the signal.

Now, all of this is fine and dandy but it doesn't  work power signals which are signals that have infinite energy but finite (average) power, that is,
$$\mathcal P_T = \frac{1}{2T}\int_{-T}^T |x(t)|^2\, \mathrm dt, \tag{3}$$ which will be recognized as the average power delivered by $x(t)$ (into a $1\Omega$ resistor) over the $2T$-second interval $[-T,T]$, approaches a finite limit $\bar{\mathcal P}$ as $T \to \infty$.  The power signals that everyone is familiar with are periodic signals (generally represented by Fourier series) and to accommodate these, we include Dirac deltas or impulses into our theory and pretend that the right side of Eq. $(1)$ "converges" to an impulse  or a sum of impulses.

  The Fourier transform of a periodic signal with fundamental frequency $f_0$ consists of impulses at the harmonics $nf_0, n \in \mathbb Z$, of the fundamental frequency.  Specifically, if the periodic signal $x(t)$ has Fourier series $\displaystyle\sum_{n=-\infty}^\infty c_n \exp(j2\pi nf_0 t)$, then the Fourier transform $X(f)$ of $x(t)$ is defined to be
  $$X(f) = \sum_{n=-\infty}^\infty c_n \delta(f-nf_0),$$
  where $\delta(\cdot)$ denotes an impulse or Dirac delta.

This has become second nature to us by now and we even incorporate impulses into tables of Fourier transforms and the like, and get into arguments about whether the Fourier transform of an impulse train is another impulse train or not. Note that the power spectral density of $x(t)$ is defined to be $S_x(f) = |X(f)|^2$ and for the case of periodic $x(t)$, we have that
$$x(t) = \sum_{n=-\infty}^\infty c_n \exp(j2\pi nf_0 t)
\implies X(f) = \sum_{n=-\infty}^\infty c_n \delta(f-nf_0)\\ ~\text{and}~ \\|X(f)|^2 = \sum_{n=-\infty}^\infty |c_n|^2 \delta(f-nf_0).$$
If you choose to "square" the sum $X(f)$ to arrive at $|X(f)|^2$, please remember that $\delta(f-nf_0)\delta(f-mf_0)$ is $\delta(f-nf_0)$ if $n$ equals $m$ and $0$ if $n\neq m$.
The autocorrelation function of $x(t)$ is the inverse Fourier transform of the power spectral density and is thus
$$r_x(t) = \sum_{n=-\infty}^\infty |c_n|^2 \exp(j2\pi nf_0 t)$$ which is also a periodic function.

But what if $x(t)$ is a power signal but is not a periodic signal and so we can't brute-force use Eq. $(1)$ and hand-wave our way to a power spectral density consisting of impulses only?  Well, one way to proceed is to re-consider Eq. $(3)$ and note that $\mathcal P_T$ can be expressed as
$$\mathcal P_T = \frac{1}{2T}\int_{-T}^T |x_T(t)|^2\, \mathrm dt = \frac{1}{2T}\int_{-\infty}^\infty |x_T(t)|^2\, \mathrm dt\tag{4}$$
where $x_T(t)$ defined as
$$x_T(t) = \begin{cases}x(t), & -T \leq t \leq T,\\0, &\text{otherwise,}\end{cases}.\tag{5}$$
Note that $x_T(t)$ is a finite-energy signal no matter how large $T$ is (as long as $T$ is finite) and thus it has a Fourier transform $X_T(f)$. Note also that $x_T(t)$ has finite support $[-T,T]$ while the support of $X_T(f)$ is the entire frequency axis $-\infty < f < \infty$. Next, note that  Parseval's relation lets us re-write Eq.$(4)$ as
$$\mathcal P_T = \frac{1}{2T}\int_{-\infty}^\infty |X_T(f)|^2\, \mathrm df\tag{6}$$
leading to
\begin{align}
\bar{\mathcal P} &= \lim_{T \to \infty}\mathcal P_T\\
&= \lim_{T \to \infty} \frac{1}{2T}\int_{-\infty}^\infty |X_T(f)|^2\, \mathrm df\\
&= \int_{-\infty}^\infty \lim_{T \to \infty} \frac{1}{2T}|X_T(f)|^2\, \mathrm df &{\scriptstyle{\text{since}~\bar{\mathcal P}~\text{is finite by assumption.}}}\tag{7}
\end{align}
But the average power is just the area under the power spectral density curve, that is, $$\bar{\mathcal P} =
\int_{-\infty}^\infty S_x(f)\, \mathrm df$$ and thus the power spectral density of this finite-power signal is defined as
\begin{align}S_x(f) &= \lim_{T \to \infty} \frac{1}{2T}|X_T(f)|^2\tag{8}\\
&= \lim_{T \to \infty} \frac{1}{2T}\left|\int_{-\infty}^\infty x_T(t)\exp(-j2\pi ft) \,\mathrm dt\right|^2\tag{9}\\
&= \lim_{T \to \infty} \frac{1}{2T}\left|\int_{-T}^T x(t)\exp(-j2\pi ft) \,\mathrm dt\right|^2\tag{10}
\end{align}
which (except for the expectation operator -- not needed because everything is deterministic here) looks a lot like the second definition of the power spectral density in the OP's question.

The autocorrelation function of $x(t)$, which is the inverse Fourier transform of the power spectral density $S_x(f)$ in Eq. $(8)$ is
\begin{align}r_x(t) &= \lim_{T \to \infty} \frac{1}{2T}\int_{-\infty}^\infty x_T(\tau)x_T(t+\tau) \,\mathrm d\tau\\
&= \lim_{T \to \infty} \frac{1}{2T}\int_{-T}^T x_T(\tau)x_T(t+\tau) \,\mathrm d\tau
&{\scriptstyle{\text{since}~x_T(\tau)=0~\text{when }~|\tau|>T}}\\
&= \lim_{T \to \infty} \frac{1}{2T}r_{x_T}(t).\tag{11}
\end{align}
Note that $r_{x_T}(t)$ has support $[-2T,2T]$ but is a poor approximation to $r_x(t)$ when $|t| > T$ because the overlap between the support $[-T, T]$ of $x_T(\tau)$ and the support of $x_T(t+\tau)$ is small.

OK, but what about random processes which is what the OP was asking about? Well, the problem is our model for  a random process is a collection of random variables $\{\mathscr X(t)\colon t\in \mathbb R\}$  whereas what we observe with our spectrum analyzers and oscilloscopes is a finite segment, say, $x_T(t)$ of one sample path $x(t)$ among the many sample paths of the process.  So, for each $t \in [-T,T]$, we know the value that the random variable $\mathscr X(t)$ took on for the one outcome $\omega$ in the underlying sample space $\Omega$ but it is quite difficult to get much information about even the  basic properties of $x_T(t)$ (e.g. is $x_T(t)$ a continuous function? does it have a Fourier transform? etc.) from the bald description of the random process as "a collection of random variables". So, let's assume that the random process is wide-sense-stationary and that it satisfies various ergodic  theorems that allow for estimation of the power spectral density and the autocorrelation function of the process from the segment $x_T(t)$ (and its Fourier transform $X_T(f)$) which is the only part of the single sample path $x(t)$ available to us. The Wiener-Khinchin theorem says the power spectral density $S_{\mathscr X}(f)$ of the process $\{\mathscr{X}(t)\}$ is given by
\begin{align}
S_{\mathscr X}(f) &= E\left[\lim_{T \to \infty} \frac{1}{2T}\left|\int_{-T}^T \mathscr{X}(t)\exp(-j2\pi ft) \,\mathrm dt\right|^2\right]\tag{12}
\end{align}
and we can approximate this by
$$\frac{1}{2T}\left|\int_{-T}^T x(t)\exp(-j2\pi ft) \,\mathrm dt\right|^2
= \frac{1}{2T}\left|\int_{-T}^T x_T(t)\exp(-j2\pi ft) \,\mathrm dt\right|^2.$$
It would appear reasonable that we should be able to express the power spectral density $S_{\mathscr X}(f)$ as the Fourier transform of the process autocorrelation function $R_{\mathscr X}(t)$ which we can approximate by $r_{x_T}(t)$??:

$$S_{\mathscr X}(f) \approx \frac{1}{2T}\int_{-\infty}^\infty r_{x_T}(t)\exp(-j2\pi ft) \, \mathrm dt??
\tag{13}$$
Well, $r_{x_T}(t)$ is known only for $|t| \leq 2T$ and might not be all that great an estimate of the value of $R_{\mathscr X}(t)$ when $|t| \geq T$. So, the Wiener-Khinchin theorem also says that the


finding the fundamental frequency in frequency domain I want to get the fundamental frequency of a signal. I used a time domain approach first. It just sums up the differences between the signal (lets say 2048 samples) and the delayed version of the same signal. I think it's not quite what autocorrelation does because I don't pad the signal with zeros or do circular convolution. Instead, I just accessed older samples. So when compairing the signals with 1 sample delay, I check samples 1 to 2048, for a 5 sample delay I check samples 5 to 2053 and so on. For the full 2048 sample delay, I check 2048 to 4096.
This works our really well because I just need the delay where the difference is the lowest. That might not be the fundamental but haveing a closer match is more important for my actual usecase. (writing an oscilloscope)

I ran into performance issues, so I wanted to use FFT instead of the time domain approach. However, I've noticed that because of the zero padding of the signa, I get some kind of triangle window. I guess what would be the equivalent in the time domain approach is using zeros instead of older samples in the delayed signal. (I've attached a screen from Matlab)

My question is if there is some way to use FFT for processing but to get rid of this "triangle window" effect. Or in other words: Is there a way to implement my time domain approach using FFT?

Cheers

A triangle is the convolution of a rectangle against itself.  Zero padding creates a rectangular window of data.

To remove the triangle (or other windowing) artifact, you can compensate the autocorrelation of your data by using the inverse of this window overlap versus lag function.


Low frequency (1Hz) High pass filter I am trying to get a useable wave from a pulse oximeter, but there is a pesky DC bias that keeps changing every now and then, so I cannot simply subtract the DC bias from it and get useable results for long. I also tried to take an appropriate sampled running average, and subtract that from my wave, and it gave me much better results than before.
However, I am told that a DSP high pass filter will work even better, and I am trying to implement that. I am not an electrical engineer and do not have much knowledge about filters, and absolutely no knowledge of DSP. I found this great website that lets me design a DSP filter and gives me a C++ library to use that filter in my code.
The problem is, I want to design a filter that cuts below 1Hz and lets through any frequencies above 1Hz or 1.3Hz max. My oximeter wave is 1.4Hz. But the design tool I mentioned above seems to indicate that it is not possible.
Any help in this would be highly appreciated.
[Edit]
I have finally managed to make it work a bit, but the DC bias is still there. (red line is the data filtered from the filter from that website, blue line is moving average subtracted from incoming signal)

And this is the frequency-gain graph of the filter

but there is a pesky DC bias that keeps changing every now and then

Fix that, if you can.  If it changes suddenly, then any filter you have is going to have odd results after the change, until it settles again.
So this is the algorithm you want (note the pseudo-C pseudo code -- data types are up to you):
static state_t state = 0;
static const cutoff_frequency = 1.0;
static const state_t gain = cutoff_frequency / (2 * PI * sampling_rate);

input_t highpass(input_t input) {
    input_t retval = input - state;
    state += gain * retval;
    return retval;
}

state_t is up to you -- it's easiest to choose 64-bit or longer floating point; excepting issues with speed or space that's almost certain to work.  Anything else and you'll need to worry about truncation.  I assume that input_t is some integer type -- if it's unsigned, you'll need to do some casting so the first line of the function works.
What this is doing is making a highpass filter by making a 1st-order infinite-impulse response lowpass filter of the input data, with unity gain, and subtracting it from the input.  As long as you have enough precision in your variables, the second line of the function insures that the output really is servoed to an average value of zero.
It's up to you to make sure that nothing overflows for extreme values of the input, and that state_t and gain fit together so that the second line of the function does not underflow (i.e., if the input changes by its smallest increment, state still "sees" the difference in retval after multiplied by gain.


Noise cancellation for known sources Creating a noise cancelling earphones that cancel/suppress speech is a huge challenge , maybe even impossible.

But what happens if we change the goals a bit:

We need to cancel the noise(music,speech) coming from a tv in the same room.
We can receive what audio the tv will create , some time before it's created through the speaker.

Does it make it possible to create great noise cancellation ? how would it work ?

If you have audio that will be played by TV beforehand, then this problem becomes almost similar to echo cancellation, i.e. adaptive filters could be used and cancellation will be great.

What we have is the audio that will be played through TV, this audio will be used as reference for adaptive filter and input to the adaptive filter would be the capture(noisy music in earphone/headphone).
so adaptive filter will try to estimate the TV's audio in the earphone capture, then this estimate could be subtracted from earphone capture to get noise free(TV noise) music in earphone/headphone.


Subcarrier spacing in OFDM The idea of OFDM is to spread data across multiple overlapping subcarriers. As they are orthogonal it is still possible to demodulate them. That means, all other subcarriers are 0 at the evaluation point of a specific subcarrier.

Orthogonality is achieved by generating a rectangular subcarrier pulse. By applying a IFFT a $\sin(x)/x$ emerges which allows orthogonality.

Another requirement for orthogonality concerns the subcarrier spacing. It must be equal to the reciprocal of the symbol period. Howver, I do not understand why. I thought it is most important to create these $\sin(x)/x$ that are shifted across the x axis on an unique position for each subcarrier.

Question:
Why would orthogonality break if the subcarrier spacing was not reciprocal of the symbol period?

There are many intepretation of OFDM to answer your question. One way I find intuitive is looking at the resolution of (Discrete) Fourier Transform.

By capturing a signal during time $T$, and by assuming periodicity outside $T$, DTFT is able to provide a frequency resolution $1/T$. Thus if you choose a subcarrier spacing $\Delta f$, you will want a frequency resolution at least $\Delta f$. As a consequence, signal must be captured at least $T = 1/\Delta f$.

One may argue that I capture then $T > 1/\Delta f$. It is fine, but the extra capturing time may catch the next symbol and cause ISI. Assuming that in a good system design, e.g. with CP or zero guard time, so that ISI is avoided, the resolution of DTFT (and DFT/IFFT) $1/T < \Delta f$. The issue here is that if $\Delta f$ is not a multiple of $1/T$, the output of IFFT will be at frequencies $k/T, k \in \mathbb{Z}$ outside the subcarrier positions and this output will be sinc-interpolation of data.

So, at the end the requirement is $\Delta f=k/T, k \in \mathbb{N_+}$ or $T=k/\Delta f$. For the sake of economy, we choose the smallest $T$ possible at $k=1$ yeilding the famous thumb rule "subcarrier spacing is reciprocal of the symbol period".


Why multi-path channel has linear phase within the coherence bandwidth? I read some sentence that said "It is usually assumed that the multi-path channel has a constant gain and linear phase within the coherence bandwidth".
I can't understand what it mean " linear phase ".
To my best of knowledge, difference path delay would caused difference phase shift in the receiving signal, but it's would be rarely that different path delay is linear, or it should be said that there is almost no such situation. Therefore, I am very curious about what it mean about "linear phase".
Can someone please explain it?
Thanks in advanced.

For the classical baseband channel model for multipath, $h(t)=\sum_ka_k\delta(t-\tau_k)$, the frequency response $H(f)=\sum_k a_ke^{-j2\pi f\tau_k}$. If $K=1$, single tap the phase would have been linear $\phi(f)=2\pi f\tau_k$. But for different path incident on the receiver, there is no guarantee that phase is linear. But let us see how phase changes between first path $k$ and last path $p$. For these incident paths, $2\pi f\tau_k-2\pi f\tau_p \approx \pi$ would cause fading. So if your frequency $f$ changes from $f_0$ to $f_0 +\frac{1}{2(\tau_k-\tau_p)}$, you will encounter fading. Hence $f_w=\frac{1}{2(\tau_k-\tau_p)}$ is called Coherence Bandwidth. It is the amount of bandwidth for which frequency response is 'flat' or constant. If the delay spread $\Delta \tau = (\tau_k-\tau_p)$ was less, $H(f)$ would be varying less over frequencies. If the difference in delay of paths was much much smaller than symbol time, that is, $\Delta \tau << T$, the width of the channel would be much much greater than width of symbols. That is $\frac{1}{\Delta\tau} >> \frac{1}{T}$. So your transmission data would see a 'flat' or constant channel in terms of frequency response. A constant channel in frequency corresponds to a single tap as explained earlier. It's phase response is linear $\phi(f)=2\pi f\tau_k$.


Rayleigh BER in Simulink I'm using the Rayleigh Channel in Simulink but I can't calculate the BER. I always get values near to 0.5 which doesn't make sense (even using high EbNo values). [In this simulink model][1] I found a trick to calculate the BER but I also would like to use delays and attenuation vectors.

Any idea of how to calculate the BER? I have tried also with a Matlab script comparing bit to bit but I also get values near 0.5.

Edit1:
I can calculate the BER but it is always around 0.5, regardless EbNo values. I'm using this simulink model.

Edit2:
I'm using delay and attenuation vector so I have problems with the array dimensions:

After the AWGN block you need to divide with the channel coefficient value before you apply the BPSK demodulator.

If you click on the Rayleigh block it gives the option to enable a "Complex Gain Path Port". You can use this output to accomplish the division mentioned above.

Instead of division you can also multiply the output of the AWGN block with the complex-conjugate of the channel gain.


What is an 'Oriented' Band-pass Filter? In this research paper, Section 4.1 talks about a filter called "Oriented Band-pass Filter".

What does it actually mean by 'Orientation'?

Also, the article gives a function,

$$H(u, v) = \frac{1}{1 + 0.414(\sqrt{\frac{u^*}{D_h} + \frac{v^*}{D_v} })^{2n}}$$

Can you please explain this function in layman's terms?

What are $u$, $v$?

I think that by orientation they mean rotation of coordinates.
$u,v$ are the horizontal and vertical frequencies while $u^*,v^*$ are the rotated frequencies.
If i understand correctly, they are looking for a line in an image with unknown rotation (aka slope) so they create of filter-bank of oriented filter, each with a different rotation.


3D point cloud to 2.5D graycale image (metric) [Halcon?] I have a problem ragarding transferring a point cloud to a 2.5D metric grayscale image (gray value gives depth information). I start with a three channel image of 32bit floats with each pixel representing the three coordinates in mm (metric).

Halcon is the preferred library, but not necessary.

I hope that someone can help me.

Greetings from Stuttgart, Germany

Max

Your 3 channels already give you the 3D coordinates, and from that 3 images, you already know which pixels locations the 3D data are projected. So you don't need to do anything. If you like to get a depth image, just set:
$$
depth(i,j) = z(i,j)
$$
In Halcon, you just select the Z-channel using access_channel.
If you desire a distance image, then:
$$
distance(i,j) = \lVert (x,y,z)^T \rVert
$$
The norm defines the distance metric you would like to use. Traditionally, I would say Euclidean $\|\cdot\|=\sqrt{x^2+y^2+z^2}$. You could do this in any language by using operations like addition, multiplication and square roots.


Lower peak to average power ratio by phase shifting A particular, frequency-domain sequence of +1 and -1  as follows:

$X=[-1, -1, +1, +1, +1, +1, +1, +1, +1, -1, +1, -1, -1, +1, +1, -1, -1, -1, +1, +1, -1, -1, -1, -1, +1, -1, +1, -1, +1, -1, -1, +1]$

is being transformed to the time-domain and transmitted over a wireless channel. It was stated that if $X$ is phase shifted i.e. multiplied by the following sequence

$w= [1, j,-1,-j,1, j,-1,-j,1, j,-1,-j,1, j,-1,-j,1, j,-1,-j,1, j,-1,-j..]$

The resulting time-domain sequence would have a lower peak to average power ratio.

Can anyone tell me why?

Each symbol (in the narrow sense - bit) in your frame modulates a sine wave of frequency depends on symbol's number and phase depends of symbol's value. It can be watched from the Inverse Fourier transform. Then these sin waves are added to each other to create time domain signal. If you change phase relations in your data frame, you change phase relations in sine waves in the sum. By this multiplication you actually do this. So obviously it alter PAPR. OFDM precoding technique (e.g. based on Hadamard orthogonal transform) is an example of using this Fourier property.


How are phase values able to capture motion from video? I know that the phase spectrum contains most of the structural information about the image
importance of phase. But I want to know more about importance of phase spectrum related to video signals.

I have read that temporal variations of phase values are able to capture most of the dynamical characteristics of the video sequence like global motion in the video phase based modelling, but I don't understand how it does?

Please consider any real time example say video of rotating wheel or wave (or video of traffic on road) if I compute its phase spectrum using Fourier Transform, the phase values captures motion of rotating wheel (or complex motion of moving car on road), but I don't understand how it does? By which property of Fourier transform could anybody can explain it to me?

Also can anybody explain the mathematical relation between motion and phase? Please correct  if I am getting wrong.

If you move an impulse in the spatial domain (perhaps representing a pixel on some object in a scene), the unwrapped phase of all the spectra associated with that point will be rotated in the frequency domain by an amount proportional to the distance of the move and the frequency.


Problems filtering image with prewitt I'm trying to filter an image vertically to get the edges, in MATLAB, but I get very different results from convultion and correlation.

CODE:

vf = fspecial('prewitt');

iyconv = imfilter(im, vf, 'replicate', 'conv');
iycorr = imfilter(im, vf, 'replicate', 'corr');

result in the image below,

why does the convultion only outputs the lower part of the mouth, and on contrary the correlation teh upper?

Prewitt kernel is rotated 180 degree between correlation and convolution.

You can check this link: http://www.equalis.com/blogpost/731635/129540/Convolution-and-Correlation-in-Image-Processing


Characterizing an unknown LTI system What are the best approaches to characterize an unknown discrete time LTI system? I believe one of the approaches is to input a known input and measure the output to find a transfer function. What are some more approaches to better characterize this black box system?

One option from the realms of adaptive filtering is the Least Mean Squares (LMS) filter depicted below:

The idea is you take the output of the unknown system, compare it with the output of your adaptive filter and minimize the difference by tweaking the filter coefficients, using a LMS algorithm. When the error $e(n)$ is zero (or more often, lower than a threshold you define) you found your unknown system's response.

There are countless resources where you can learn about LMS. In this case, given

the input vector: $$X(i) = [x(i) \space x(i-1) \ldots x(i-N+1)]^T$$ and the coefficients vector: $$W = [w_0 \space w_1 \ldots w_{N-1}]^T$$

The output of the adaptive filter is: $$\hat{y}(i)=X(i)^TW$$

Now define the cost function to be: $$e(i)^2 = \big(d(n)-\hat{y}(i)\big)^2$$
It can be shown that the gradient of the cost function with respect to the coefficients is:$$\hat{\nabla}(i) = \frac{\partial e(i)^2}{\partial W} = -2e(i)X(i)$$
Since the gradient points towards the maximum of the cost function you update the coefficients in the opposite direction: $$\hat{W}(i+1) = \hat{W}(i)-\mu\hat{\nabla}(i) = \hat{W}(i)+2\mu e(i)X(i)$$
where $\mu$ controls the adaptation speed. I've ignored a lot of details but this is the idea.

Note that the filter will produce a frequency response that is adapted to you input signal $x(n)$. It's not guaranteed for any other frequency. For example, in the response depicted below, my $x(n)$ is composed of 3 tones. You can see only at those particular frequencies does the filter (N=3,6,10) match the true response (H1(z)). So consider feeding it a wideband signal like white noise.

If you want to learn more you can check for example:

Wikipedia article

B. Widrow, Adaptive Signal Processing, Prentice-Hall, 1985


Can someone explain what a butterfly is in regards to DSP? I often see the term "butterfly" in discussion and methods relating to DSP. What is a "Butterfly?"

Wikipedia has this one covered:

  In the context of fast Fourier transform algorithms, a butterfly is a portion of the computation that combines the results of smaller discrete Fourier transforms (DFTs) into a larger DFT, or vice versa (breaking a larger DFT up into subtransforms). The name "butterfly" comes from the shape of the data-flow diagram in the radix-2 case, as described below.1 The earliest occurrence in print of the term is thought to be in a 1969 MIT technical report.[2][3] The same structure can also be found in the Viterbi algorithm, used for finding the most likely sequence of hidden states.

(Source)


information filter instead of kalman filter approach I read many sources about kalman filter, yet no about the other approach to filtering, where canonical parametrization instead of moments parametrization is used.
So I would like to learn on examples how to use this filter with information matrix and information vector (not mean and covariance from KF).
Can anyone help me with this, give some examples, webpages, or parts of implementation code to be a basis to learn?
Really appreciate!

Information filters are an application of a learning algorithm. One example of an information filter would be a recommendation system. On the other hand, a kalman filter is restricted to linear tracking. Its not learning per se.

If you are looking for machine learning/AI based approaches to SLAM, then you should probably look at Particle Filters.


Lower frequencies stronger in FFT? I'm calculating FFT from microphone input. I notice that lower frequencies always seems to have more power (higher dB) than higher frequencies.

I cut the data into frames of 24576 bytes (4096*6).
Apply Hamming window: input[i] *= (0.54d - 0.46d*(double) Math.Cos((2d*Math.PI*i)/fs));
Run it through FFTW Process1D().
Convert from complex numbers: output[i] = 10.0 * Math.Log10((fout[i * 2] * fout[i * 2]) + (fout[i * 2 + 1] * fout[i * 2 + 1]));
Average out 6 values to get a complete FFT of 4096 bytes.
Paint pretty picture (colors mapped to a palette).

The averaging (pt. 5) is done to lower FFT noise.

As the image shows both with sound and with mic off there is more energy (and more noise) in the lower frequencies. This indicates its more than just a mic/sub problem.

My questions:
 1. Is this expected? Why?
 2. Any standard way of fixing this? Almost looks like some Math.Tan() magic could lift it up where it belongs.

My goal is to be able to identify the top 5 frequencies without lower frequencies winning by default.

Yes, this is very much expected. What you see is a "pink" spectrum, i.e. constant energy per relative bandwidth, as compared to "white" which is constant energy per absolute bandwidth. For pink signals the energy between 1-2kHz is the same as from 2-4kHz (each represents a doubling of bandwidth or an "octave".

Most natural audio signals (Speech, music, movies, etc.) have a pink-like spectrum. Also most acoustic background noises (microphone noise, HVAC noise, air pressure fluctuations, generic background) tend to be pink or even more skewed towards lower frequencies.

The human auditory systems works that way, too. In the inner ear the audio signals are broken down into bands of constant relative bandwidth (called "Critical Bands") which are roughly the same as a third octave spectrum.

The best way to look at audio spectrum data is to plot it on a logarithmic frequency scale.


What are the advantages of barker sequence over marker sequnce? I read that barker sequence is used for synchronization among other uses.  Barker sequences are rare and they have special properties - this part I understood.  I am missing some thing when I think about how they score over simple marker sequences for synchronization.  You want to know where the information bits start - start with a marker sequence preferably which does not occur in the information.  And you can synchronize.  How do barker sequences score over marker sequences?

A big advantage of Barker Codes over a generic marker is strong correlation when the codes are aligned and very low correlation for all other shifts, even by one sample.  This offers increased resistance to bit errors that would otherwise cause a false detection.
The presence of a Barker sequence can be optimally detected using correlation, which is the process of multiplying bit by bit the received sequence to the reference Barker code sequence at a given alignment and summing the result of the multiplications (multiply and accumulate). With a Barker sequence, you will get maximum correlation when the received sequence and reference sequence are aligned, and near 0 for all other shifts. This is an ideal characteristic for synchronization, and not a property of generic marker sequences.
Depending on the noise resistance desired, longer codes can be chosen, and/or code sequences can be repeated. Like the use of a generic marker sequence, the Barker code itself would not be an allowed pattern in the data payload, or properly escaped.
As an example, 10110111000 is an 11 chip Barker Code. As an easy way to explain the correlation as a multiply and accumulate, replace all the 0's in the sequence with -1, and all 1's in the sequence with +1. Therefore when correlating the multiplications before summing are an x-nor: 0x0 =1, 0x1=0, 1x0=0, 1x1=0. For demonstration, consider a sequence of three repeating Barker codes for a header:
received sequence
101101110001011011100010110111000
As the sequence is received, we will correlate it by shifting it against our reference sequence and performing the multiply and accumulate for each shift.
Case 1 aligned:
When a reference code is aligned with the received code, the correlation (multiply and sum) is:
101101110001011...
10110111000
1x1=+1
0x0=+1
1x1=+1
1x1=+1
0x0=+1
1x1=+1
1x1=+1
1x1=+1
0x0=+1
0x0=+1
0x0=+1
Sum=+11

Case 2 misaligned by 1 sample:
01011011100010...
10110111000
1x0=-1
0x1=-1
1x0=-1
1x1=+1
0x1=-1
1x0=-1
1x1=+1
1x1=+1
0x1=-1
0x0=+1
0x0=+1
Sum = -1

The same result of -1 occurs for all shifts, except when the code aligns at which point the correlation is maximum (the length of the code).  Since the length of the code is odd, the correlation cannot sum to the ideal value of 0, but gets quite close! Also note that the code can be inverted, and in that case it would sum to -11 when aligned (for a maximum distance of 22 between positive and negative copies of the code).
Maximum Length Sequences (generated with Linear Feedback Shift Registers) also have this shift-correlation property and include significantly longer sequences than Barker codes, which finds application in GPS and Direct Sequence Spread Spectrum communications where signals are received below the noise floor and recovered using these correlation techniques.
Another view from coding theory is the concept of Hamming distance.
(see [https://en.wikipedia.org/wiki/Hamming_distance] 1) To be consistent for that approach we would need to repeat what I did above with symbols 0 and 1 and use an x-or to do the multiply. If the symbols in two sequences are different a 1 is added, if they are the same 0 is added (x-or). The Hamming distance is the total addition of the x-or results between the two sequences.  If every symbol didn't match (maximum negative correlation) the Hamming distance would be 11, and if every symbol matched, the Hamming distance would be 0 (for the 11 chip Barker sequence I used).  Therefore the Hamming distance between 01001000111 and 10110111000 is 11 (maximum distance possible), and the Hamming distance between the code and cyclically shifted versions of the code would be approximately half way in between: For example, the Hamming distance between 01001000111 and 10010001110 = 5. In fact, since the Barker code has noise-like properties, the likelihood is highest for the distance to be near the middle for all random patterns such as what may be represented by our data payload.  It is much more likely that the Barker code will have distances between 4 and 8 to patterns in our data but much less so from 0-3 and 9-11; but that still can an will happen, especially with such a short sequence, unless we encode our data to ensure those sequences do not occur unless bit errors existed. The concept of Hamming distance, and choosing code words of maximum distance is a foundation of error detection and correction.  The bigger the distance, the more errors we can correct (see the reference for specific details on all of that).


MPEG-7 low-level audio descriptors I'm working on an application which need to use the full set of 17 MPEG-7 low-level descriptors. Is there any library I can use for it? Preferably with a Python wrapper/API...

I just did some research on existing implementations of the MPEG-7 audio descriptors and found the following implementations:

The official implementation of the 17 low-level descriptors (part 1 and part 2) (MATLAB).

Alghough the descriptors are implemented in MATLAB, it should not be a big issue to implement them in Python, based on the MATLAB code. Maybe the following resources are also of interest:

The Timbre Toolbox, which implements a number of features described in this JASA article, which include some of the MPEG-7 audio descriptors (MATLAB) and
An an online computation of all 17 audio descriptors with the MPEG-7 Low Level Descriptor Extractor of the TU Berlin.


In image processing, what is the difference or relationship between noise and artifact? I would like to know the theoretical difference or relationship between noise and artifact in images, I barely know the relation between noise and distributions, but I'd like a deep explanation if it is possible related to medical imaging field.

Do you have a link, paper at hand?

The difference is that noise may obscure features in an image, while artefacts appear to be features but are not.  If the 'problem' is structured, it is probably an artefact, whereas if it is random, it is probably noise (as a generalisation).

A computed tomography example:  noise will make the image look grainy, and make small differences in contrast difficult (or impossible) to identify.  A streak artefact on the other hand has structure to it, and looks like the patient has a region of low density where they actually don't.  Likewise with ring artefacts - they look like the patient has ring shaped structures within them.  Artefacts are generally caused by assumptions made in the development of the reconstruction algorithm (though not always).

It is similar to the difference between noise and interference in some other fields - interference is structured in some way, having some particular source, while noise is less predictable.

EDIT: Just found something in the glossary of "Computed Tomography" by Willi A. Kalender, 2011 (978-3-89578-317-3):

artifact: part of the contents of an image that does not have a counterpart in the physical object being imaged; in CT images an artifact can appear as false structures (e.g. aliasing artifact, beam hardening artifact, motion artifact, partial volume artifact) or as a corruption of the CT numbers.
noise: contributions to a measured signal due to random processes; in principle these contributions are of statistical nature and do not carry any information about the signal.

If that's of any use to you.


Effects of interchanging sine terms with cosine terms Suppose we have a real signal $x(t)$. Now, we know that $x(t)$ can be represented as a sum of sines and cosines.
w be the angular frequency.

If $a(\omega)$ be the coefficients of the cosine terms, and $b(\omega)$ be the coefficients of the sine terms, then if we perform a transformation on this signal such that now $b(\omega)$ are the coefficients of not sines but cosines and similarly, $a(\omega)$ the coefficients of the sines, then how will the waveform transform$[x(t)]$  look like?

What will be the similarity between this transformed waveform and the original waveform?
What is the name for this transformation (is it Hilbert transformation? )

Please attach graphs for an example $x(t)$ to clarify.

If you have a $T$-periodic signal $x(t)$ with a Fourier series expansion

$$x(t)=a_0+\sum_{n=1}^{\infty}\left[a_n\cos(n\omega_0t)+b_n\sin(n\omega_0t)\right],\quad \omega_0=\frac{2\pi}{T}\tag{1}$$

then its Hilbert transform is given by

$$\mathcal{H}\{x(t)\}=\hat{x}(t)=\sum_{n=1}^{\infty}\left[a_n\sin(n\omega_0t)-b_n\cos(n\omega_0t)\right]\tag{2}$$

because $\mathcal{H}\{\cos(t)\}=\sin(t)$, $\mathcal{H}\{\sin(t)\}=-\cos(t)$, and $\mathcal{H}\{c\}=0$ for constant $c$.

This is not exactly the same as just swapping the coefficients $a_n$ and $b_n$. The signal with swapped coefficients (assuming $a_0=0$) is

$$y(t)=\sum_{n=1}^{\infty}\left[a_n\sin(n\omega_0t)+b_n\cos(n\omega_0t)\right]$$

which is related to $\hat{x}(t)$ by

$$y(t)=-\hat{x}(-t)\tag{3}$$

It is also equal to the Hilbert transform of $x(-t)$:

$$y(t)=\mathcal{H}\{x(-t)\}\tag{4}$$

So swapping the coefficients is equivalent to applying a linear time-varying system to the signal. It cannot be achieved by a Hilbert transform alone (which is a time-invariant operation).

Note that for even signals ($b_n=0$) and for odd signals ($a_n=0$) one does indeed obtain the Hilbert transform (for odd signals with a negative sign) by swapping coefficients. This is shown in the example in hops's answer.


Is spectral leakage due to windowing 'different' for the DTFT and DFT? I'm currently trying to improve my foundational understanding of spectral leakage due to windowing. The more I have been reading about it, the more I'm starting to think that 'windowing spectral leakage' can mean two completely different things depending on whether you're talking about the DTFT or DFT. I'm hoping to get some confirmation/clarification on what I discuss.

With regards to the DTFT, this transform is essentially the Fourier transform equivalent for a discrete signal. This transform assumes that the signal is periodic with that period being infinity. With regards to windowing, spectral leakage is induced due to the fact that you are multiplying a discrete signal $x[n]$, by some window function $w[n]$. The corresponding convolution in the frequency domain causes the spectrum of $X[f]$ to be smeared or 'leaked'.
Now moving onto the DFT, the idea of spectral leakage induced by windowing seems to be a completely different concept. Unlike the DTFT, the DFT assumes the signal is periodic (with that period being less than infinity). In this context, it is often stated that spectral leakage will occur if the collected signal does not contain an integer multiple number of points to the digital frequencies contained (i.e. if you have a $3\textrm{ Hz}$ signal with sample rate $10$ samples/sec, and collect $18$ samples, you will have spectral leakage). Thus if you recorded a sinusoid at interval where it's not able to 'finish' a complete cycle, the next period will restart and there will be a discontinuity in the signal. While this makes sense that you would be inducing 'false' frequencies or leaking frequencies, it seems to be a completely different idea from the convolving concept of the DTFT.

My primary question:

Is 'window spectral leakage' completely dependent on context (DTFT vs DFT)?
Furthermore, is it safe to say that the spectral leakage induced by each case separate from one another? As in, spectral leakage from the DTFT is induced by frequency convolution, and spectral leakage from the DFT induced by the assumption of periodicity? Or is there something that ties these two things together?

The concept of spectral leakage is not dependent on the context, and it's the same thing for the DTFT as it is for the DFT. It's probably helpful to remember that the DFT of a finite length sequence equals a sampled version of the DTFT of the same sequence:

$$\begin{align}X(\omega)&=\sum_{n=0}^{N-1}x[n]e^{-jn\omega}\qquad\text{(DTFT)}\\
\tilde{X}[k]&=\sum_{n=0}^{N-1}x[n]e^{-jn\frac{2\pi k}{N}}\qquad\text{(DFT)}\end{align}$$

Obviously we have

$$\tilde{X}[k]=X\left(\frac{2\pi k}{N}\right)$$

Now assume we have two sinusoidal signals windowed with a rectangular window such that in one case there's an integer number of periods inside the window, and in the other cases there isn't:

In both cases we get spectral leakage due to the windowing process. However, in the first case, the samples of the DTFT computed by the DFT are all at the zeros of the DTFT (except for one sample at the sinusoid's frequency). This is not the case for the second signal. That means that in the first case the DFT doesn't "see" the spectral leakage due to the sampling process, but in the other case the spectral leakage is also visible in the DFT.

These two cases are shown in the figure below. The DTFT is shown in blue, and the corresponding DFT is in red:


Stochastic approximation algorithm The goal is to find the FIR filter coefficients $\mathbf{h} = [5;3]$ with the help of the adaptive FIR filter $\mathbf{w}$ of order $p = 2$.

I have implemented the Stochastic approximation algorithm in MATLAB for this adaptive system identification problem.

$\mathbf{w}[k] = \mathbf{w}[k-1] + (1/k) \cdot \mathbf{x}[k] e[k]$

As input sequence I use white gaussian noise with zero mean and variance of 1.
Measurement noise is also WGN with variance of 1.

$N = 1000 $ ... number of measured datapoints

For the input sequence above my algorithm works fine and converges and I find the right $\mathbf{w}$.

However if I change my input sequence to $x[k] = 0.1 \cdot \cos(\pi k /(N/10))$ my algorithm doesn't find the right coefficients.

So my question is:

Is the Stochastic approximation algorithm restricted to special input sequences, and if so what is the optimal input sequence to use?

The issue is possibly that the input signal you have chosen is not persistently exciting.  This means that the signal doesn't "excite" enough modes of the filter in order to be able to accurately estimate its parameters.  Another way to think about it is that it doesn't have enough energy in enough places in the spectrum: just at the frequency of the cosine, whereas white noise has energy everywhere in the spectrum.

From the paper:

  Recent work (e.g. [1–5]) has shown the convergence and robustness properties of adaptive estimation
  and control algorithms to be closely related to persistence of excitation of regression vectors (equivalently
  that they have a sufficient number of spectral lines [5])


Mapping Saw Wave That Is Between -1 And 1 To A Sine Wave I am trying to map values of a saw-wave that goes from -1 minimum value to 1 maximum value as a periodic signal onto a sine-wave.
Repeating the matrix of this waveform that is 257 element in length oscillates as expected as a saw-wave, but if I change the function into a sine-wave -1 maps to 0, then 1 also maps to 0, and when two zeros are together when we loop the waveform with "repmat" matlab function, there appears discontinuities in the sine-wave function...How can I map this saw-wave, so the sine-wave perfectly oscillates with all the values used from the saw-wave that is between -1 and 1 values?

f=(-1:1/128:1)';

c1=repmat(f,1000,1);

c2=sin(pi*c1);   % -1 of saw wave maps to 0, 1 maps also to 0

sound(c2,44100); % two adjacent 0's cause discontinuities

I also tried the fft method which gives me the needed coefficients but if my waveform is complex this would become more complicated to adapt:

f=(-1:1/128:1)';

m=fft(f);

m(3:256)=0;
m(1)=0;

m2=ifft(m);

m3=m2/max(abs(m2));

m4=repmat(m3,1000,1);

sound(m4,44100);

You can scale your sawtooth before you put it into the sine function, like so:

f=(-1:1/128:1)';

ramp = f * (length(f)-1)/length(f);  % Scale down to not include the last 1 into the sine
c1=repmat(ramp,1000,1);

c2=sin(pi*c1);   % -1 of saw wave maps to 0, 1 maps also to 0

plot([c1 c2]);
%sound(c2, 44100);


Downsampling a high pass filter If I have an ideal high pass filter as the one in the picture:

And I downsample it by 2, e.g. $h_{new}[n] = h[2n]$, is my understanding correct that the resulting filter would look like this in the frequency domain:

Your understanding in not correct now. And it can be seen from the equation which states the relation between DTFTs of $h[n]$ and $h_{new}[n] = h[2n]$ as follows :
$$ H_{new}(e^{j\omega}) = \frac{1}{2} \left( H(e^{j\frac{\omega}{2}}) + H(e^{j\frac{\omega-2\pi}{2}}) \right) $$ The general form of this relation for $h_{new} = h[Mn]$ is

$$H_{new}(e^{j\omega}) = \frac{1}{M} \sum_{k=0}^{M-1} H(e^{j\frac{\omega+2\pi k}{M}}) $$

It's up to you to sketch the resulting expanded and shifted spectrums that add up to create the final spectrum, to see that it's a lowpass filter with cutoff frequency of $\omega_c = \pi/2$

The following MATLAB/Octave code also helps to see it practically:

h = fir1(256, 0.75, 'high');    % Create a HIHPASS filter with cutoff wc = 3*pi/4
h2 = h(1:2:end);                % Create h_new[n] = h[2n] , downsampled filter

figure, plot(linspace(-1,1,1024), abs(fftshift(fft(h,1024))));
hold on
plot(linspace(-1,1,1024), abs(fftshift(fft(h2,1024))),'r');
title('DTFT magnitudes for h[n] and h2[n] = h[2n]');
legend('h[n]','h[2n]')

The resulting plot is:


Eigenvalues of a vector I'm trying to calculate (with Matlab) the minimal and maximal Eigenvalues $\lambda_{min}$ and $\lambda_{max}$ of the autocorrelationmatrix $R$ to adjust the stepsize of an adaptive filter.

My input signal $x$ is a vector of colored noise.
And my plan is to do something like

lmin = min(eig(xcorr(x)));
lmax = max(eig(xcorr(x)));

But it fails becaus 'For standard eigenproblem EIG(A), A must be square.'

I think this is a little bit strange, because in my opinion the output of the autocorrelation is always just a vector.

Please help me to find the point, whether I need to use something else than Matlabs 'eig' function or if it's mathematically not possible to calculate the eigenvalues of a single vector.

The autocorrelation matrix is a Toeplitz symmetric matrix that can be formed as:

Rxx = toeplitz(xcorr(x));

and that should give you a sensible (square) matrix that can be passed to eig.


sampling a continuous time signal using MATLAB? I have a continuous time signal $x(t)=cos(2 πft)$ where $0\leq t \leq 2$.

I want to plot discrete time signal $x(n)$ where $0 \leq n \leq 19$ for sampling frequency 100 Hz and signal frequency 10 Hz.

My MATLAB code is follow:

clc
clear all
close all
t=0:0.1:2
f=10
x=cos(2*pi*f*t)
subplot(211)
plot(t,x)
n=0:19
fs=100
Ts=0.01
x_a=cos(2*pi*f*n*0.01)
subplot(212)
stem(n,x_a)``

Please let me know, am I correct?

This information can be found other places as well but I will step through it here using MATLAB.

You say your sampling rate is $100$ Hz and the signal is $10$ Hz. Lets define those along with the sampling period ($1/100$ seconds).

sampleRate = 100;
samplePeriod = 1/sampleRate;
signalFreq = 10;

Now you want $2$ seconds of this, so define your sampled time vector.

nT = 0:samplePeriod:2;

Now we create the signal, and plot the result.

signal = cos(2*pi*signalFreq*nT);

figure
stem(nT, signal)
grid on
xlabel('Time (sec)')
ylabel('Amplitude')

If you only want to plot $0 \leq n \leq 19$, then you may instead set nT = [0:19]*samplePeriod;.


Properties of Stochastic Processes I have a basic question about stochastic processes:

When some informations such as wss, uncorraleted sampled, white about random signal (say x[n]) are given, what do we exactly have?

For example white information:
I know if x[n] is white, its power density function $ P_{x}(e^{j\omega}) $ is constant value independent of $ \omega $. I suppose it means autocorrelation function of x[n] consists of dirac? (I am not sure here, If I am wrong, please correct me).
Uncorraleted samples: I know if the random signal is uncorrelated, its covariance equals to zero which means its autocorrelation equals to square of its mean. But uncorraleted samples means this, or something else?
What about wss?: I know if the signal is wss, mean (average) of x[n] is constant which is useful information when $ x_{c}(t) $ (constinuous time form of x[n]) and its means are given and we can say directly mean of x[n] is the same as mean of $ x_{c}(t) $. Moreover, we have $ \phi_{x}[m]  =  \phi_{x_{c}}[mT] $ where T is the sampling period and $ \mu_{x_{c}}^{2} = \mu_{x}^{2} $ where $ \mu $ is the power of the signals (second moment).

Any correction, addition, comment would be appreciated.

Pass through an example may be more convenient:

Suppose we have a system like this:

where, x[n]: wss, $ \eta_{x} = 0 $ (mean value), uncorrelated samples, $ \sigma_{x}^{2} $ and h[n] are given.

How to calculate $ \eta_{y}[m] $ the mean, $ \phi_{y}[m] $ the autocorrelation and $ \mu_{y}[m] $ the average power of y[n]?
How to calculate $ \eta_{v}[m] $ the mean, $ \phi_{v}[m] $ the autocorrelation and $ \mu_{v}[m] $ the average power of v[n]?

As pointed out in the comments, this question is rather broad, so here I will just answer your last question about the autocorrelation of the output signal $y[n]$. You might know the following relation concerning the power spectra of input and output sequences:

$$S_y(\omega)=S_x(\omega)|H(\omega)|^2\tag{1}$$

where $H(\omega)$ is the frequency response of the LTI system. From the inverse Fourier transform of (1) you get a corresponding relation between the autocorrelation functions:

$$\phi_y[m]=\phi_x[m]\star h[m]\star h^*[-m]\tag{2}$$

where $h[m]$ is the system's impulse response, and $\star$ denotes convolution.


What's the best way to compare the same object from different "crops" or angles within multiple images? I read a similar question, but I'm not sure the techniques discussed will help me.

Take for example the two following images:

I'm trying to determine if the above images are sourced from the same image as well as whether the objects in the following image are the same or similar enough to safely consider them the same.

Thanks in advance
Joe

There are many ways to achieve this. You could use an image based, direct similarity metric such as Sum of Absolute Distances, Sum of Squared Distances or Normalized Cross Correlation. For a general overview look here. These approaches generalize to template matching where a convolutional procedure is conducted, such as the one here.

Template matching methods could tolerate rigid transformations (rotation, translation) and scaling. However, for the case of affine deformations (e.g. the second image set in the question), they would not suffice. It is then better to benefit from feature based matching techniques, where some repeatable keypoints are first extracted from the image, and invariant descriptors of the local neighborhoods are computed and store. This is done for both images. Finally the descriptor are compared using a nearest neighbor search. SIFT is one of those robust descriptors.

Note that, while in the first scenario, the template transformation can be defined simply by a translation, the second case is not like that. You should either solve for affine parameters, or compute a general perspective homography to represent the pose of the second image w.r.t the first one.

The problem is very general, and often times analyzed under the roof of object detection. You could review this literature and choose a method suited to your problems. OpenCV implements most of the techniques, I described in this answer.


Are there characteristics of photographic film that would make it an analog medium? There was an article a few years ago entitled "Analog is Not the Opposite of Digital".

Film is not analog, period. I used to shoot film on a Canon AE-1 from the ’70s. Now I have a digital SLR from Canon, and they’re obviously extremely different. But we have to be careful not to confuse ‘old’ and ‘new’, with two very specific terms like analog and digital.
The word digital, to most people, refers to a device that can capture, store, or display data in a binary fashion. Ones and zeros, on and off, digital is all about numbers. Digital shouldn’t be confused with binary, of course, as digital simply means concrete values. Any system that utilizes solid values (or digits) is digital, binary is simply the most common system. Digital cameras capture light with a sensor, that light is converted into data (numbers), so the use of the word ‘digital’ for your cell phone camera or DSLR is accurate.
Analog, however, is a very abused word. I would venture a guess that the significant amount of readers have used the word ‘analog’ to refer to film cameras. If the new, fancy robot cameras are ‘digital’ then our aging film cameras are ‘analog’, right?  Not at all. Older cameras capture light with film, which is basically plastic, gelatin, and silver halide. When you take a photo (perhaps of your dog drinking a beer), photons hit this material and produce a latent (invisible) image, that can later be brought into view by bathing the film in various chemicals. You could write hundreds of blog posts on film development alone, but the point is that film photography is a chemical process.

I understand the author's point about this misuse of the word "analog" in general, but it seems as though chemical film possesses many of the characteristics of an analog medium, in that it can capture a continuous spectrum of color values over a particular range.
Does this characteristic (or any other) of chemical film make it directly comparable to other analog media such as a record or a cassette, or does the fact that the method of storage of the light information is not a continuous wave like a phonograph record preclude it from being classified as such?
Since I know it will likely show up in a search for these topics, I had found Analog Photography appears to be a subfield of the art of photography, in which "progressively changing [the] recording medium" creates the image, regardless of whether chemical film or digital capture is used.  The technique itself seems to be the analog process in that case, and not any sort of statement about the photographic medium.

Film isn't absolutely "analog", as in continuous.  Every individual silver halide molecule, after exposure and development, is either metalized or not; and there are a finite number of these molecules in every frame of film, thus quantizing the exposure measurement.  However the density and location of the film grains and silver halide molecules is semi-randomized, which helps noise shape or dither the sampling, quantization and aliasing noise more than if they were in a regular grid with fixed step values.

So some current differences are not that film is in some sense a truly continuous measurement system (at least down to the Planck level), but more might be that the "sampling" is dithered, and also that the dynamic range may be higher than in the most common digital formats (8 bit integer per channel).