texts: What is the LOW frequency resolution rule analogous to Nyquist?

texts:
If I'm analyzing a signal in the frequency domain, I know about the well known Nyquist criteria that the sample frequency must be > 2x of the highest component present in the signal.

However, there are very few references related to the lowest frequency I can analyze.    This problem comes up in a lot of natural data where we have limited windows of observation because the instruments to do so were only invented or in wide spread use in the last 50-100 years, and natural cycles can be decades, centuries, or longer.

One would think that you could analyze a component where at least one period was in the sample window, but I find because of required windowing to prevent edge effects that the signal often disappears (it depends on the phase).  I then thought it was a mirror of Nyquist and I could resolve any signal with at least 2 periods in a window, but actual analysis shows that to not be true when there's a another overlying signal whose period is similar.

I found through experiment that I needed at least 5 periods to resolve a signal where there's another signal near the same frequency.

Is there a formal way of representing this problem or a definitive reference I can use?

Would the same issue apply to wavelet analysis in the same manner?


texts: Can FFT convolution be faster than direct convolution for signals of large sizes?

texts:
Let's say I have a 1D signal of size $N$ and am trying to filter it with a 10-tap FIR filter mask. When $N$ is large, the number of multiply-accumulates would approximately equal
$$2 \times 10 \times N = 20N\quad\text{floating point operations}$$
This convolution algorithm is thus $\mathcal O(N)$ since the computation scales linearly with the input data size.
However, when would it make sense to use an FFT convolution to find this convolution? If an FFT is $\mathcal O(N\log N)$, wouldn't this mean FFT convolution does not scale as well as direct convolution? If this is the case, how can FFT convolution be faster than direct convolution?


texts: LowPass Filter & Window *Img & Code* - Frequency Response has no lobe roll off

texts:
Not sure if I am doing something silly but I am not seeing lobe roll off on the frequency response of my windows and the low pass filter I made in MATLAB.
Please can you help ... several hours just looking at this and my textbook -.-

It should be like this showing lobe roll off

%% Low Pass Filter
Fs=10*10^3; cutoff = 1000; %100 Hz
fc = cutoff / Fs ;
b = 800/Fs ;%Transition bandwidth of 800 Hz
M = 51 ; %Length of Filter
n = 0:1:M-1;
blackman_window = 0.42 - 0.50 * cos((2*pi*n)/(M-1)) + 0.08*cos( (4*pi*n) / (M-1));
hamming_window  = 0.54 - 0.46 * cos((2*pi.*n)/(M-1));
hanning_window  = 0.50 - 0.50 * cos((2*pi.*n)/(M-1));
blackman_window_ft = fftshift(fft(blackman_window)/M);
hamming_window_ft  = fftshift(fft(hamming_window)/M);
hanning_window_ft  = fftshift(fft(hanning_window)/M);

sinc_func = sin(2*pi*fc.*(n-((M-1)/2))) ./ (2*pi*fc.* (n-((M-1)/2)));
sinc_func((M+1)/2) = 1 ;
sinc_func_ft = fftshift(fft(sinc_func)/M);

h_lpf = sinc_func .* blackman_window;
h_lpf = h_lpf ./ sum(h_lpf); %Normalise
H_lpf = fft(h_lpf);
H_lpf = fftshift(H_lpf);

figure(1);cla;clf
subplot(3,2,1);hold on;
plot(blackman_window,'DisplayName','Blackman')
plot(hamming_window,'DisplayName','Hamming')
plot(hanning_window,'DisplayName','Hanning')
plot(sinc_func,'DisplayName','Sinc')
legend
title('Windows and Sinc Impulse Response');hold off

subplot(3,2,2);hold on;
plot(20*log10(abs(blackman_window_ft)),'DisplayName','Blackman')
plot(20*log10(abs(hamming_window_ft)),'DisplayName','Hamming')
plot(20*log10(abs(hanning_window_ft)),'DisplayName','Hanning')
plot(20*log10(abs(sinc_func_ft )),'DisplayName','Sinc')
legend
title('Window and Sinc Frequency Response (dB)');hold off;

subplot(3,2,3);hold on;
plot(h_lpf);
title('Low Pass Filter Impulse Response');hold off;

subplot(3,2,5);hold on;
plot(abs(H_lpf),'DisplayName','Magnitude');
plot(real(H_lpf),'DisplayName','Real');
plot(imag(H_lpf),'DisplayName','Imag');
legend
title('Low Pass Filter Frequency Response');hold off;

subplot(3,2,6);
plot(20*log10(abs(H_lpf)));
title('Low Pass Filter Frequency Response (dB)');


texts: Confusion regarding discrete signal value?

texts: I know that discrete signals have finite precision. But does that means that those signals can have only integer values(for example,2,3,4 etc) but not decimal values(for example 2.1,3.2,4.2 etc)?


texts: Approaches to Segmentation in images of water surface

texts:
this is my first question here. Therefore please help me inform me when I miss any good practice in asking questions here.

Given a short video sequence of an water surface, I would like to detect if there are any objects visible. The video sequence is taken from an moving air vehicle.

An example frame looks like this :

Currently I only process one frame.
Assumption: most of the image is water
I found that water can look very different, even in the same frame.
Approach A: Looking for edges with the sobel-operator. Con: waves create edges, too. Therefor correct thresholds are difficult to find.
Approach B: For every pixel measure distance to greater local neighbourhood (for example 128x128). Con: Finding correct thresholds for scene
Note: I used different colourspaces (rgb,yuv,hsv). I found hsv to be usefull for approach B. Hue of water seems to be nearly constant. Value of most objects is quite different from water (even waves). Saturation varies with waves and is not suited for thresholding.

I would like to ask:

What are approaches that might work finding non-water-regions/pixel?
Any thought to my approaches? Are there usal improvements that I might have missed?
Is it probably that combining the results of neighbouring frames improves detection or reduces false-positiv-detections?

Thank you for your inputs.

Edit:

here is another image that I find difficult.

It has no relevant colour information. Therefor I only present the lightness channel.


texts: What are the differences between estimation and detection optimal filters (matched and Wiener)?

texts: What is the differences between the two of them? I know that the matched filter is a detection filter because it detect what is the peak of the signal and return the time delay of reflected and transmitted signal. Whereas the Wiener filter is an estimation filter because it estimates the mean-square-error estimation of the degraded signal. But I'm not too sure on what is the definition of the detection filter itself, and the estimation filter definition. Any help would be appreciated.


texts: How to establish transfer function of a speaker?

texts:
I'm working on a signal processing project about vocal system, and I'm trying to use controlling theories to solve the problem.

I need to get the transfer function, $H(s)$, of a speaker, from electrical signal to voice. But it seems that I have to use a microphone to pick up the voice. Considering the microphone itself has a transfer function $G(s)$, I will finally get $G(s)H(s)$ rather than $H(s)$.

So, is there any method to get TF $H(s)$ from electrical signal to vocal signal?

Well I have another problem. I've seen someone else estimating TF $G(s)H(s)$ by experiment. The bode diagram looks so strange and complex that I think it will be difficult to get a mathematical description of the system. So what method is usually used in controller-design when the object has a complex freq. response?


texts: how to detect and capture scene change in a recorded video

texts:
I have a project to detect and capture scene change in a recorded video.Actually its a recorded video of a godown (of about 2 weeks) and I have to detect and capture pic whenever someone comes in.

I have tried vlc player motion detector. It detects changes very well but does not capture the detected scenes.
Is there any way that can help?

thanks and regards


texts: Can I create my IF signal on SDR before bringing data to PC?

texts:
I'm relatively new to SDR's and I'm trying to determine what sort of BW I need for my laptop to get the best resolution possible for an FMCW system. I have an X310 Ettus with a UBX160 daughtercard which should give about 1 meter of resolution theoretically. Anywhere close to this will probably work just fine for my application. But at the moment I only have a PC capable of a 1 Gig ethernet hookup. The fact is though this will be a pretty short range system, and I expect a pretty narrow range of results (approx bandwidth of a few hundred kHz). Is there a way to mix the echo with the FMCW waveform on board the ettus and just deliver the sub MHz band to the laptop? Or do I really need to shell out all that cash to get the range resolution this unit could deliver?
Sorry this is a newb question, I've been playing around with GNU Radio and uhd a little and I'm leaning towards expecting a "no" on this question - but would love some input from the pros out there. Thanks in advance!


texts:
I want to ask Question about the Fourier series in continuous time domain.
I am following signal and systems 2nd Edition by Alan Oppenheim.
I have confusion in understanding the statement that
Specifically, suppose that $x(t)$ is real and can be represent in the form 3.25. then since $x^*(t) = x(t)$, we obtain

$$x(t) = \sum^{+\infty}_{k\ =\ -\infty} a^*_k e^{-jk\omega_0t}$$

Then it means the equation 3.25 is for both Real and imaginary?
Equation 3.25

$$x(t) = \sum^{+\infty}_{k\ =\ -\infty} a_k e^{jk\omega_0t}$$


texts:
I was just wondering... The formula I learned to calculate the energy of the signal is expressed in the time domain:

$$E_x^{\text{time}} = \sum_{n=-\infty}^{\infty} |x[n]|^2$$

Then, what does the amount of energy gotten from the magnitude spectrum mean?

$$E_x^{\text{frequency}} = \sum_{f=0}^{fs/2} |X[f]|^2$$

Supposing that both amplitudes are expressed in decibels, $fs/2$ is the maximum representable frequency (by Nyquist theorem), and $E_x^{\text{frequency}}$ is taken over the same time frame as $E_x^{\text{time}}$. Should the quantities be the same?


texts: SVD Based Noise Filter on the Samples in Toeplitz Form

texts:
Given a sequence of length $2N$
$$
x[1], x[2], \dots, x[N], \dots x[2N]
$$
We can construct the following matrix along with its SVD.
$$
M =
\begin{pmatrix}
x[1] & x[2] & \cdots & x[N] \\
x[2] & x[3] & \cdots & x[N+1] \\
\vdots  & \vdots  & \ddots & \vdots  \\
x[N+1] & x[N+2] & \cdots & x[2N]
\end{pmatrix} = U S V^T
$$
Small singular values can be removed and matrix recomputed. Input signal can be replaced with e.g. 1st and last cols of $M$.
Unfortunately, I do not know the name and motivation behind this procedure.
Where can I find a reference with a detailed description of such a filter?
What is matrix $M$ called?
Example:

Test code in Mathematica:
(* define test signal and noise *)
length = 2^10 ;
amplitude = 1.0 ;
signal = amplitude*Sin[2*Pi*0.2741*Range[2*length]] ;
SeedRandom[1] ;
level = 0.1 ;
noise = RandomVariate[NormalDistribution[0.0,level*amplitude],2*length] ;
input = signal + noise ;
(* filter *)
matrix = Most[Partition[input, length, 1]];
{u,s,v} =SingularValueDecomposition[matrix,Tolerance -> 10^-15 ] ;
keep = 2 ;
s = Take[s,keep,keep] ;
u = Take[u,All,keep] ;
v= Take[v,All, keep] ;
matrix = u.s.Transpose[v] ;
result = {First[matrix],Rest[Last[matrix]]} //Flatten;
(* compare *)
Grid[{
{
ListPlot[signal,ImageSize ->300],
ListPlot[input,ImageSize ->300],
ListPlot[result,ImageSize ->300]
},
{
Periodogram[signal,Automatic,Automatic,KaiserWindow[#,Pi/2]&,PlotRange -> All, ImageSize ->300],
Periodogram[input,Automatic,Automatic,KaiserWindow[#,Pi/2]&,PlotRange -> All, ImageSize ->300],
Periodogram[result,Automatic,Automatic,KaiserWindow[#,Pi/2]&,PlotRange -> All, ImageSize ->300]
}
}]


texts: Is there a window with both a lower side lobe level and higher side lobe roll of rate than Blackman

texts: Is there a window with both a lower side lobe level and higher side lobe roll of rate than the Blackman window - even if it is poor or useless for some other reason?


texts: Continuous-time mathematical formula for deconvolution filters

texts:
I have an elementary function $p:\mathbb{R}^2\rightarrow\mathbb{R}$ which (locally) represents an image. It's a polynomial, and its the result of the following 2D convolution:

$$p=f\star G\star \mathbf{1}_{A}$$

where $G(x,y)$ is a 2D Gaussian, and $\mathbf{1}_A$ 'box' function:

$$G(x,y)=e^{-\frac{(x^2+y^2)}{2\sigma^2}}, \quad\mathbf{1}_A(x,y)=1 \text{ for }x \text{ and } y < a, 0\text{ otherwise.}$$

Is there a closed form for a function with which I can convolve $p$, which results in $f$?

Currently I am 'avoiding' the use of such a filter by performing Fourier deconvolution, but mathematically this opens a huge can of worms.


texts: Algorithm for Discrete Signal Decomposition

texts:
I am building a sensor and I am trying to understand how to process the signal that it generates. The sensor has a library of reference signals. When 'event a' occurs, it produces signal A, when 'event B' occurs it produces signal B. EX:

It is possible for multiple (but few) events to occur, and their signals would add together linearly. I am looking for an algorithm for signal decomposition which takes advantage of the fact that the decomposition is discrete (It can give integer multiples of the library signals only). I also need some measure of the reliability of the identification.

Ex:

  Input:

     (I made this image by taking A + 2 * B + noise)

Example Output:
A + 2 B  -  95% confidence
A + B    -  3% confidence
A        -  1% confidence
None of the Above - 1% confidence

I've noticed that I can compare two different possibilities by taking an inner product of the normalized proposed solution and input signals. For instance:

Input Signal .  (A + 2 B) = 0.999

but

Input Signal .  (A + B) = 0.985

So the first solution is a closer match to the input signal. I thought about using a Newton-Raphson to maximize this inner product, but when the actual sensor is in use it could have 1000s of signals in it's library and that would not be realistic.


texts: Numerically Controlled Oscillator (NCO) for phasor implementation?

texts:
If I want to simulate a digital oscillator (or phasor) to modulate an arbitrarily length signal:
$$y(t) = \cos(2\pi f_ct)\quad\text{where}\quad t = \frac{n}{f_s}\quad\text{for sample}\quad n$$

What is the advantage of using a numerically controlled oscillator (Answered here) versus simply incrementing $n$ over the length of the input?

main():
  n = 0                 # total samples processed
  N = len(input)        # num samples in input
  blksz = N / 16        # don't create phasor all in one call
  # calculate phasor over multiple passes
  while ( n < N )
      update_phasor(blksz)
# calculate phasor
update_phasor(blksz):
  for (i = 0; i < blksz; i++)
      y[i] = exp(j*2*pi*fc*n/fs)
      n++

The NCO referenced above and the pseudocode implementations give different results and I am trying to understand which makes more sense.

Also, what is the best way to synthesize the phasor if the desired frequency needs to change on the fly?


texts:
I have a signal $s(t)$ of duration $T$ with limited bandwidth. In matlab, I use a sampling frequency $fs$ and I have $L=T*fs$. After, I have $r=ones(size(s))+s$. To have a limited bandwidth, I filter with a raised cosine filter: $rr=filter(h,1,r)$ with $h$ the raised cosine filter.

My question: What is the duration for the signal $rr$? Personnally I would have said $T$. How to prove that?


texts: Image processing, recognizing a small feature in a larger image

texts:
I am trying to write an image processing program to recognize bubbles in oil.
It has been suggested I try computing the convolution of the image and an image of a typical bubble.
i.e. : ifft(fft(a).fft(b)) will have high peaks near bubbles.

where:
ifft : inverse fourier transform
fft : fourier transform
a : the image
b : kernel image (bubbles)

Is this a good way to solve this problem? Also, for the dot product to work, I need the images to be the same dimension. This is not likely to be the case.


texts:
I am confused trying to understand the Proof of Fourier Transform from Oppenheim book Signals and Systems. I am pasting the equations directly from the book:

  $$\widetilde{x}(t)=\sum_{k=-\infty}^{+\infty}a_ke^{jk\omega_0t} \tag{4.3} $$
  $$a_k=\frac{1}{T}\int_{-T/2}^{T/2}\widetilde{x}(t)e^{-jk\omega_0t}dt \tag{4.4} $$

  where $\omega_0=2\pi/T$. Since $\widetilde{x}(t)=x(t)$ for $|t| < T/2$, and also, since $x(t)=0$ outside this interval, eq. $(4.4)$ can be
  rewritten as

  $$a_k=\frac{1}{T}\int_{-T/2}^{T/2}x(t)e^{-jk\omega_0t}dt=\frac{1}{T}\int_{-\infty}^{+\infty}x(t)e^{-jk\omega_0t}dt$$

My question is In the last equation on the L.H.S, there is no
$\widetilde{x}(t)$ but it was mentioned in the eq. $(4.4)$
And then we have: $$Ta_k= \left.\frac{2 \sin(\omega T_1)}{\omega}\right\rvert_{\omega =
k\omega_0}$$ That is, with $\omega$ thought of as a continuous variable, the function ($2\sin\omega T_1/ \omega$) represents envelope of $Ta_k$.
Here I want to know what is meant by envelope?I want to add more details from the book exactly " Specifically we have plot $T_0a_k$ rather than $a_k$ and we have also modified the horizontal spacing in each plot. here is the plot "

c) I want to know that how $T_0a_k$ will be drawn in continuous?


texts:
I am new to EE area and I met a question in my homework:

Can someone help me a little bit? I think it is quite easy!

Thanks!!


texts: Removing the effects of camera motion in a video

texts: I have frames of a video consisting of both a moving camera(translation, rotation, zooming in and out) as well as moving objects. I need to recover only the trajectory of objects which means I need to align the images spatially by removing camera motion. Can anyone guide me on how I can do that?


texts: Find highest frequency of a very manky signal

texts:
I'm trying to compute the highest frequency (as can be sampled) in some pretty manky looking discrete time-dependent signals. My current method - a discrete fourier analysis - fails for some pretty awful looking but clearly oscillating signals (with discernable highest frequencies).

My current method is to compute the discrete fourier transform, locate the local maxima, and perform weighted averages around these peaks (to undo the 'smearing' of the dft across the discrete sample frequency bins).
Here it is in MATLAB.

% compute spectrum
dft = abs(fft(weight));         % real response
dft = dft(1:floor(end/2) + 1);  % keep only pos freqs
dft = dft/sum(dft);             % normalise

% find present modes modes
mode = zeros(0, 2);

[pks, locs, widths, ~] = findpeaks(dft);
for k=1:length(locs)

    % only consider significant peaks
    if pks(k) > 0.1

        av_width = min([floor(peak_width_factor*widths(k)), locs(k)-1, length(dft)-locs(k)]);
        inds     = (locs(k)-av_width):(locs(k)+av_width);
        av_ind   = (inds * dft(inds)) / sum(dft(inds));
        av_freq  = (av_ind - 1)/T;

        mode = [mode;
                av_freq, pks(k)];
    end
end

This produces a matrix of present modes (of a spectral significance above 0.1/1) where each row is the frequency and significance of the mode.

This works great for signals like this (time signal left, dft right with the frequencies of detected modes labeled):

but fails for signals like these (where we see and expect frequencies close to 2; not an order of magnitude smaller!)

Often I expect a certain frequency (e.g f=2.6 in the above signals) and can judge myself the mode is present by the average period between the local maxima. I tried codifying this - computing the average time between local maxima - in MATLAB, but it was pretty unreliable:

[~, peak_inds] = findpeaks(weight);
cycle_periods = diff(peak_inds) * dt;
av_cycle_period = sum(cycle_periods)/length(cycle_periods);
av_freq = 1/av_cycle_period;

My signals are well enough sampled (around 15 to 20 values per observable manky period) to resolve these manky modes visually.
I've studied time series and random processes at an undergraduate level, but we never really went too deep into spectral analysis. So:

How can I reliably compute these highest-frequency modes in my very 'unfourier' signals?
Why does my current DFT analysis incorrectly deduce very low frequency modes in these manky signals?
Why would, if multiple present modes have dissimilar frequencies, my naive distance-between-local-maxima-average method fail at extracting the highest frequency?


texts: Filtering a signal to have a constant exponential moving average

texts:
Problem

I wish to filter a discrete-time input signal, $x[n]$, so as the output signal's, $y[n]$, exponential moving average, $a[n+1] = a[n] + \alpha(y[n] - a[n])$, tracks around a constant reference $\mu_x$ such that $\mathbb{E}\{y[n]\} \approx \mu_x$, and also an approximately constant variance, $\text{var}\{y[n]\} \approx \sigma_x^2$, while maintaining as much information about the incoming signal as possible.

Naive Motivating Solution

1) Maintain an average of the input, $b[n+1] = b[n] + \alpha(x[n] - b[n])$ and an approximation of the signal variance, $c[n]$.

2) Define the output to be $y[n] = \frac{x[n] - b[n]}{\sqrt{c[n]}}$ which tracks a reference of 0 with unity variance.

Question

There are a number of problems with this approach which make it not very numerically robust, such as when low values of variance occur.

Is there a more intelligent way of doing this? Perhaps using a dynamic linear filter? Or high pass filter?

Other Points

The inputs can be considered independent, and are drawn from one of several dynamic distributions.

The application, if you are interested, is comparing the amplitudes of recent inputs with each other by assigning them a "reward" output. This is high reward for relatively higher input values and a low for lower lower values. However, I desire the mean reward to be an approximately constant output mean. Furthermore, there may be some medium-time dynamics in the inputs, such as decaying toward zero over time.


texts:
I have particle activity as shown in the left pane of animation below. The activity is clustered and it moves slowly. Sometimes these clusters merges together. On the right side of it, its shuffled version is shown. Each frame of the activity on the left is shuffled by a fixed permutation matrix P.

The problem I have is the inverse problem. Often I record what looks like activity on the right. Shuffling is allowed. My hypothesis is that there is a permutation $P^{-1}$ (fixed for each frame) which can transform activity on the right to the left.

Not sure if this has a unique solution or the problem is even well-posed. I am not familiar with literature on this. Any tip, hints or reference would be of great help.

PS: Are there some features which remains invariant post permutation?


texts:
Given the following poles and zeros map I have to identify qualitatively, which filter type is represented (Low-Pass, High-Pass, ...). In the solution it says, that this is a high pass filter.

However, I plotted the amplitude response (see below) and to me this looks more like a notch filter or a band stop.

Where is my mistake in reasoning?

Thank You!


texts:
I am a bit confused about the Kelly Lochbaum Ladder Filter. Image below taken from here.

My question is, is this an LTI (linear time-invariant) filter? I thought that it was LTI, but I am not sure anymore.


texts: Standard deviation is varying with change in amplitude of DC signal

texts:
I am trying to make sense of the below code. I am taking a large sample so that I get good estimates for standard deviation.

N = 25;
h1 = ones(1, N)./N;
L = 1000000;
s = ones(1, L);
n = 0.1*randn(1, L);
x = s + n;
stdx = std(x);
y = filter(h1,1,x);
stdy = std(y);

This code will give me stdx = 0.1 and stdy = 0.02, which is expected since $\sigma_y=\sigma_x/{\sqrt{N}}$as we are essentially doing an average filter.

But giving more amplitude to the signal  s = 10*ones(1, L); and running the same code above is giving me a stdy = 0.0344. Shouldn't the standard deviation be same in both the cases since the signal s has no variance?


texts: CZT: Do Bluestein algorithm perform backward transformation?

texts:
I am in discussion with one of my colleague that he says changing sign of the exponent in Bluestein's algorithm for Chirp-Z transformation, transforms the spectrum from frequency to time domain?
I am sort of confused, since they have invented ICZT very recently and it works in a different way.


texts: Barcode localization in difficult environments

texts:
I'm creating a program for localizing types of barcodes.
(I know plenty exist, that's not the question :) )
I'm focusing for the time being on 1D barcodes.

The difficulty I have is that the barcodes I need to be able to locate could be situated anywhere within an image with other information (text/images/...), inequal lighting conditions (-> gradients), and barcodes could be scaled and orientated in any way.

[Edit] Unlike most barcode-related cases however, my application does not need to be run on a smartphone, so can allow myself some processing bandwidth

I have been researching as much as possible in books and publications, but think it's time to ask the experts...

What is the most appropriate approach for reliably localising 1D barcodes in difficult environments?

Currently, I'm applying otsu thresholding, basic edge detection and hough transforms.
Problems i currently face:
- otsu thresholding is global and doesn't handle local gradients on large images well. are there localized algorithms that exist?
- hough locates the lines of the barcodes, but also other items (in particular text), and it's difficult to filter only the barcodes (thresholding the min number of pixels in the line throws out a lot, but also throws out the smaller barcodes)

Could anyone point me in the right direction?

I'm wondering about Fourier and Gabor filters, but they seem more complex, so want opinions before i start.

PS: I am doing everything in Octave and Java


texts: Output of LTI (in time and frequency $\omega$ domain) : when input goes through LPF

texts:
I would like to raise a mathematical question :
Let's say we are been given :
$$x(t) = \begin{cases} \cos(\pi t) & |t| \leq 0.5 \\
                                       0 & \textrm{otherwise} \end{cases}$$
If $x(t)$ goes through an LPF : $H(\omega) = \begin{cases} \pi - 0.5|\omega| & |\omega| \leq 1.5\pi\\
0 & \textrm{otherwise} \end{cases}$
Compute and sketch :

Fourier coefficients of $x(t)$, as well as $X(\omega)$
$y(t)$ and $Y(\omega)$ where $y$ is the output of the LTI. Is y(t) periodic? If so compute it's period and its trigonmetric fourier anlysis coefficients
I did the math and I've got this:
$a_n = \frac{2}{\pi}\frac{\cos(\frac{n\pi}{2})}{1-n^2}$ which can not be defined for $n=-1$ or $ n=1$
update:* We have to compute $a_1$ and $a_{-1}$ in seperate: If we do the computations we gain : $a_1= a_{-1}=0$
I then did :
$$c_m = \frac{1}{2}(a_m + jbm) \rightarrow c_m = \frac{a_m}{2}$$
$$x(t) = \sum_{m=-\infty}^{\infty}c_me^{jm \omega_o t} \rightarrow
X(\omega) =2\pi \sum_{m=-\infty}^{\infty}c_m\delta(\omega - m\omega_o)$$

What I tried:
Now we know for an LTI : $y(t) = x(t) * h(t) \rightarrow Y(\omega) = X(\omega)H(\omega)$
but since $H(\omega)$ exists only in $\delta=[-1.5\pi, 1.5\pi]$ and $\omega_o = \pi$ then we only need $X(\omega)$ in $\delta$ .
But , $X(\omega)$ is $0$  for $|m|=1 \rightarrow |\omega| =\pi$ and therefore $X(\omega) \neq 0 \rightarrow m=0$ in $\delta$ so :
$X(\omega) = 2\pi a_0 \delta(\omega)$ in $\delta$
we also know $\int_{-a}^{a}f(x)\delta(x)dx=f(0)$  so:
\begin{align}
y(t) &= \frac{1}{2\pi}\int_{-\infty}^{\infty}Y(\omega)e^{j\omega t}d\omega \rightarrow \int_{-1.5\pi}^{1.5\pi}\frac{1}{2\pi}Y(\omega)e^{j\omega t}d\omega \\
&\rightarrow \int_{-1.5\pi}^{1.5\pi}\frac{1}{2\pi}(\pi^2a_0 -0.5a_0|\omega|)e^{j\omega t}\delta(\omega)d\omega = \frac{1}{2\pi}\pi^2 a_0 = \frac{1}{2\pi}\pi^2 \frac{2}{\pi} = 1
\end{align}
I assume that it should have been periodic but it's not. So where did my solution go wrong?


texts: How exactly to implement the fractional part of a 'fractional equalizer'?

texts:
I have written my own LSE equalizer (that works on the bit level), and literature I have read says that 'fractional' based ones are much better. So how exactly does one get 'fractional' bits? Is it as simple as resampling the softbit decisions to 2x, 4x or whatever of their bit rate and applying the equalizer like that?

Thanks


texts: MATLAB phase of 2D rectangular pulse's Fourier transform

texts:
I'm trying to plot the graph of the phase of the Fourier transform of a 2D rectangular pulse. I've been able to evaluate the FFT but I'm not sure if the phase is correct because there are some tilts I can't explain.

What are those diagonal lines that appear in the phase?
How can I eliminate them?

Here there is the code:

clc
close all
clear
npoints=512;
perc=1;
dt=6*1E-7/(npoints*perc);  % Tempo di campionamento
df=1/(npoints*dt); % Frequenza di campionamento

t(1)=0;
f(1)=0;
for k=2:npoints/2
    t(k)=(k-1)*dt;
    t(npoints-k+2)=-t(k);
    f(k)=(k-1)*df;
    f(npoints-k+2)=-f(k);
end
t(npoints/2+1)=t(npoints/2+2)-dt;
f(npoints/2+1)=f(npoints/2+2)-df;

ts=ifftshift(t);
fs=fftshift(f);

figure
[X,Y]=meshgrid(ts,ts);

D = npoints/2;        % to indicate origin at the center of the function
a = 5;          % change it to enlarge or reduce the pulse
y = repmat(1:npoints,npoints,1);
x = y';
rect = zeros(npoints);
rect(D-a:D+a-1,D-a:D+a-1) = ones(2*a);
rect=(rect);
surf(X,Y,rect);
shading interp
axis tight
title ('Rect 3D');
rect=ifftshift(rect);
figure, surf(X,Y,rect);
shading interp
axis tight
title ('Rect 3D shifted');

R = fft2((rect));
R = fftshift(R);
[X,Y]=meshgrid(fs,fs);
figure;
surf(X,Y,abs(R));
shading interp
axis tight
title('Fourier Transform of Rectangular function');
%plot real part
figure;
surf(X,Y,real(R));
shading interp
axis tight
title('Real part');

Rm=abs(R);
imm=imag(R);
re=real(R);
re(abs(re) < 1e-8) = 0;
imm(abs(imm) < 1e-8) = 0;
R=re+imm;
phase=angle(R);

%plot phase
[X,Y]=meshgrid(fs,fs);
figure
surf(X,Y,phase);
shading flat
axis tight
title ('Phase of the rect');


texts: Should the order of a filter be proportional to the sampling rate?

texts: I have designed a fir filter with my desired frequency response at a specific sampling rate and filter order. If I were to redesign my filter for a new higher sampling rate, should I change the filter order? I am having trouble understanding how the same filter order influences the desired frequency response at different sampling rates.


texts: OFDM Preamble Ripple Structure in Frequency Domain


texts: What is the low-pass filter that gives minimal transient response?

texts:
As title says, what is the low-pass filter that gives minimal transient response? By minimal transient response, I mean minimal settling time with reference to very tiny amount of divergence of steady-state values.

I am asking for analog case, but digital case would also be fine.


texts: I was wondering if anyone knew anywhere I can obtain benchmark comparison of various algorithms used for similarity recognition. In my research I found one study from 2008 or so, but misplaced the URL.


texts: Applications or physical interpretation of auto-convolution?

texts:
I wonder if anyone has any experience with auto-convolution. In particular i'm interested in understanding the physical interpretation of it.
I understand what convolution, correlation and auto-correlation are, also i'm aware that the definition of auto-convolution will be something like
$$f\ast f = \int_{-\infty}^{\infty} f(\tau)f(t-\tau)d\tau $$
but i still don't get what are the implications o the meaning of it.  I've been looking for a while and so far i haven't found any good or detailed explanation (on constrast with auto-correlation).
So, if anyone has any experience dealing with this topic or has any intuitive interpretation that could share, i'll appreciate it. Thanks.


texts:
I want to do time-frequency analysis of a signal. So i decided to use spectrogram function from scipy. This function returns, three parameters f,t,Sxx. I am using below arguements:

x : my signal
fs : 1000 Hz
window : hann, with size = 300
nperseg : 300, same as window size
noverlap : 25

My question is about returned time segments of this function.

As far as I understand the first time segment for which we calculate fft must be 300 samples, then subtract 25 from it calculate fft for next window of 300 samples,and so on. So the signal of length 1s must be divided in 4 time segments. .275,.500,.775,1. , but the function returns :

f,t,Sxx=spectrogram(sig1,fs=1000,window=get_window('hann',300),noverlap=25,scaling='spectrum',nperseg=300)
t
Out[172]: array([ 0.15 ,  0.425,  0.7  ])

can someone please suggest if I am thinking wrong or there is a bug in spectrogram function.


texts:
I am at the moment trying to implement a contra harmonic filter,to remove the pepper noise from an image. But the implementation of it seems a bit incorrect, and I am not sure what is wrong.

void ContraHarmonic(Mat src, Mat dst, int kernel, double P)
{
  cout << "got it" << endl;
    Mat temp = src.clone();
    //copyMakeBorder(temp,temp,kernel-1,kernel-1,kernel-1,kernel-1,BORDER_CONSTANT,Scalar(0,0,0));

    for(int row = kernel/2; row < temp.rows - kernel/2-1; row++)
    {
        for(int col = kernel/2; col < temp.cols - kernel/2-1; col++)
        {

        double den=0,num=0;
          for(int i = -(kernel/2); i <= (kernel/2) ; i++)
          {
            for(int j = -(kernel/2) ; j <= (kernel/2); j++)
            {
                den += pow(temp.at<float>(row+i,col+j),P);
                num += pow(temp.at<float>(row+i,col+j),P+1);
            }
          }

        dst.at<float>(row,col) = num/den;
      }
    }
    dst.convertTo(dst,CV_8U);
}

The output i receive from it point at something being wrong with it.. I am not sure whether the implementation is incorrect, or the programming isn't done correctly.  This is the output i receive.. Which point at something being wrong.

--Updated--

Updated with the answer provided and the image looks like this


texts: FFT shows systematic deviations from analytical result

texts:
I am calculating the numerical Fourier transform of an exponential decay exp(-|t|) and compare it to the analytically calculated result, a Lorentzian. I find that the numerically calculated spectrum contains systematically larger amplitudes than the analytical one and that this deviation increases with frequency. I'm using the python library numpy.fft.
My time array contains [-100., -99.99, ..., 0.00, ..., 100.] and the signal is correspondingly
[3.72007598e-44, ..., 1.0, ..., 3.72007598e-44]. I've plotted this "original" exponential together with the back Fourier transforms below.

I then calculate the FFT (and adjust the phase to get a purely real result). When I compare the numerical result to the expected Lorentzian
$2 \cdot \frac{1}{1^2 + (2 \pi f)^2}$, I get very good agreement at low frequencies that however gets worse towards high frequencies. The numerical result (absolute value = real part) is systematically larger than the analytical one (labelled "calc" in this plot).

This deviation is of course also visible when plotting the ratio and difference of these two curves ("calc" refers to the analytic result). The ratio of the quantities seems to be independent on the step width and the maximum of the time array.

The inverse FFT of the numerical result (shown in the first plot, labelled "ift") seems to agree even a bit better with the original function than the FFT of the analytically calculated Lorentzian (labelled "icalc"). This depends however on the choice of the time array; it can also happen that the inverse FT of the Lorentzian is slightly longer close to the exponential function.
My objective is to evaluate experimental data in the frequency domain by models in the time domain. I would like to replace the usually used analytical FT (which restricts us to very idealized models) by a numerical FT. The first step is of course to show that the different FTs yield the same result... Of course, the experiment will not measure the spectrum with infinite precision, but I would very much appreciate any help to get at least 1-2 orders of magnitude better in the agreement numerical / analytical.


texts:
Internet is very much costly in India , like 3-5 USD  per GB @ 3G networks and even in broadband 100 USD for a unlimited banwith for 2 MBps but here DTH (Direct-To-Home broadcast service ) is very common and cheap , like 4$ for a month with most of the channels . I think to buffer a 480p video we would need 200 KB/s (2MBps ) plan in internet but how DTH providers give that much bandwidth for TV service for such a low price .

My main point here is :

Can we use DTH service for Internet using (direct from satellite like the DTH works).
Can we get same bandwidth for using internet just at same price just in case of DTH .


texts: How to deal with "weird" phase plots in bode diagram when designing a controller

texts:
I am trying to design a balance controller for a robot. With MATLAB simulink I arrived at the transfer function between the input and the pitch angle for the robot. I have plotted the bode and Nyquist plot below.

As you can see, the Nyquist plot looks very weird, and in the bode diagram the signals at low frequencies are attenuated. To fix this, I added a zero at the peak of the magnitude plot, so
$$\omega_i=10 \text{rad}/s \Rightarrow \tau_i = \frac{1}{\omega_i}=0.1\text{s}$$
This results in an I-term which looks like this: $C_i(s)= \frac{\tau_is+1}{\tau_is}=\frac{0.1s+1}{0.1s}$
When I add this to the controller the bode plot now becomes.

Now I also want to add a lead-term to prevent overshoots but the way I usually do that is to change the crossover frequency $\omega_c$ such that there is good phase margin, define some constant $\alpha$, and then find the lead term with the formula $\tau_d = \frac{1}{\omega_c \sqrt{\alpha}}$. Then I can make my lead-term with
$$C_d = \frac{\tau_d s + 1}{\alpha \tau_d s + 1} $$
But when I look at the current phase plot, it all around looks like there is no trouble regarding the phase margin. The phase is never close to $-180^{\circ}$ so I am not sure how to proceed with such a weird looking phase plot.
So my question is, how do I find my lead-term when the phase plot looks like this? Is there another way to do than I have proposed?
Edit
The open loop transfer function between the input and pitch angle is $G(s) = \frac{ 5.893e05 s^4 + 7.71e08 s^3 + 1.435e11 s^2 + 2.566e12 s + 7.683e09}{s^7 + 2418 s^6 + 1.317e06 s^5 + 1.977e08 s^4 + 1.312e10 s^3 + 1.796e11 s^2 - 9.658e11 s - 1.525e13}$
Using the "pole" command in MATLAB I noticed that there was a pole in the RHP which should cause some instability.


texts:
I am rather new to the world of signal processing, and am struggling to understand a fundamental concept: How are filters actually implemented?

I have read a significant portion of this online book, and scrounged the internet, finding snippets of useful information here and there, but I cannot quite tie it all together.

In brief, I have a vibration profile (~200,000 samples) in the time domain, which I would like to analyze using a particular frequency weighting. This involves a four-step filter, using a (Butterworth) high-pass and low-pass, followed by an a-v transition and an upward step. The analog equations are shown below.

Yet I am still struggling on how to implement them. I have seen mentions of the bilinear transform, used in programs in MATLAB/Python, but such implementations seem to omit the 's' or 'p' variable from the filters shown, typically creating the numerator 'B' and the denominator 'A' from the coefficients of the analog filters.

Other sources suggest taking the Fourier of the filter equations and the signal, multiplying, then applying the inverse Fourier to bring it back into the time domain. Yet I cannot discern how to take the Fourier of an equation.

My same problem occurs with using convolution in the time domain, as I have not understood how to convolve an equation yet.

I am sure I have come near the answer, and I know this is a basic principle, but I cannot seem to comprehend it, and my Googling is off the mark. Any help on this would be appreciated. General advice is also welcome, as I have a lot to learn.

** If important, I am currently working in Python, using the scipy.signal package, with access to MATLAB for testing and evaluation purposes.

Currently I am looking at substituting j2πf for p, then simply running the data through the equations and taking the inverse Fourier to bring it back into the time domain.

Alternatively,  my eyes have just noticed what was staring me in the face: the second equation given for each filter is in terms of 'f', which would suggest I am able to simply run the frequency-domain data through this equation before taking the inverse Fourier to recover the time-domain. If this is a correct realization, this question can be closed (or I can 'answer' my own, now self-evident, question).


texts: Sampling rate: camera resolution when taking a photo of an LCD

texts:
Consider a digital camera that faces an LCD and that is perfectly focused on the LCD. The LCD shows an image of $n \times n$ pixels. What resolution does the camera have to have to perfectly capture the image displayed on the LCD?

(Let's keep this simple and assume that camera pixels and LCD pixels are of the same size, pixels are "atomic" objects, no perspective distortion, etc.)

I'm asking because I have found the following statement in a paper:

  "To satisfy the Nyquist criteria for image resolution, each pixel of
  the image shown on the LCD should be sampled by 2 or more pixels in
  the camera."

I don't understand this. Intuitively, I would say that if camera and LCD were perfectly aligned, then a camera with $n \times n$ pixels would suffice. My (poor) understanding of theory also suggests the statement is wrong, because every row or column of the image on the LCD can have a maximum frequency component of $n/2$. The minimum sampling rate would be $n$, not $2n$.

What am I missing?


texts: Amplitude spectrum (transfer function) of signal?

texts: I have one question related to finding amplitude spectrum (transfer function) of signal knowing that output signal is time derivative of input signal. I have the answer graph but I don't understand the reason.


texts: Finding the Nyquist Frequency of Irregularly Sampled Discrete Data

texts:
According to L. Eyer et al. (1999), the correct Nyquist frequency of irregularly sampled data is

$$
\frac{1}{2p}
$$

Where $p$ is the largest value such that:

$$
 \forall t_i,t_i = t_1 + n_ip \\
\text{where}~n_i \in \mathbb{N}
$$

However, the paper gives no information on how to calculate $n_i$, other than stating that:

  $p$ is a kind of greatest common divisor (gcd) for all $(t_i - t_1)$.

Which doesn't really help. So, how does one calculate $p$ or $n_i$ for a given time series. Or, more importantly, is there a better way of calculating the Nyquist frequency of irregularly sampled data?

N.B. That the common methods of taking $p = s$ (the shortest sampling rate in the time series) or $p = \delta s$ (the average sampling rate in the time series) or $p = S$ (the longest sampling rate in the time series) all result in explicable behaviour (as outlined in the paper) -- mainly that super-Nyquist frequencies can be resolved (not through aliasing trickery).


texts:
I don't know if I am posting in the correct place. I wanted to ask about dilation and erosion in bitmaps. I have as a homework to implement in a programming language these morphological operations but I fail in understanding how they work. I know that erosion will shrink an object and dilation will enlarge it.

For example I have
$$A=\matrix{0 &1 &0 &1 & 1\\
0 &1 &0 &1 & 0\\
0 &1 &1 &1 & 0\\
0 &1 &1 &0 & 0\\
0 &1 &1 &1 & 0},\;
se=\matrix{1 &0 &0\\
1 &1 &0\\
0 &0 &1}$$

What would be the eroded and dilated output and how to compute it step by step?


texts: Why do these 2 methods give different solutions?

texts:
I need to solve what is underlined in red for $x_i$, nut currently I'm interested in the right side of the equation only.

On the left I sarted by doing the Laplace transform of $x_u'$ and $x_u$, and this gives the correct solution.

On the right side I tried to derive $x_u$ first, plug everything into the equation and then do the Laplace transform, and this gives a completely different solution. Is this method incorrect? Or am I doing something wrong?


texts: How do I find variance from the PSD of a stochastic process?

texts:
I have a time series that consists of noise and a signal, shown here windowed and Wiener filtered:

and the PSD of just the noise (used in filtering):

I want to find the variance of the noise using the PSD in order to estimate the significance of a measurement. I'm not super well-versed in signal processing, so I apologize in advance for any conceptual blunders! From digging around Wikipedia and StackExchange, I found that I can estimate the variance of a deterministic process like so:

$\sigma_{x} ^{2} = \int_{-\infty}^{\infty} S_{x}(f)df - \mu_{x}^{2}$

where $S_{x}(f)$ is the PSD and $\mu_{x}$ is the mean (side note: should there be a factor of $\frac{1}{\pi}$ or $\frac{1}{2\pi}$ in front of the integral?).

I then found this answer that confirmed that.

My question is: what modifications, if any, do I have to make to that formula to apply it to a stochastic process (the noise)? Also, the answer linked above specifies a stationary process- do I have to detrend the noise or can I assume that it is stationary?


texts: Bilinear transform with pre-warping for systems other than classical filters

texts:
I only seem to be able to find online information about applying bilinear transform + pre-warping to filters (like butterworth, etc.) with only one edge frequency that is purposely 'designed' into it. So the BT is easy to apply, just substitute with the usual tangent formula including edge frequency and sampling frequency.

But what if my analog rational transfer function describes a physical system (like an elastic mechanical body with several eigen-frequencies, or an electric circuit network)? Do I have to apply the BT to every pole/zero separately then? And what if the transfer function comes from measurements - I don't know the poles/zeros, because I did not 'design' the 'filter'? Does it really mean I have to factorize it (which we know is sensitive to error) just to transform it from analog to digital?

PS, Sidenote: I am going to implement the numerics involved (in C#) myself. So Matlab is certainly a nice thing, but it does not help me if not all relevant code is publicly available.


texts: Resources for modulation/demodulation techniques with example code

texts:
So I remember seeing someone recommending a DSP book that had actual code for modulation and demodulation of a communication signal. Does anyone have any recommended resources they could point me towards that might having this? I've seem lots of things that have block diagrams, but no actual implementation of how it's done. I'm interested in learning how to sync and demodulate Binary DPSK, QPSK, and chirm modulation (or sweep modulation like they use for radar and sonar I believe)

For example, I know for Binary DPSK I can use differential encoding, by taking a sample in time, and multiplying it by another sample in time from 1 bit ago. Let's say my bitrate is $R$, then my time per bit is $T_b = \frac{1}{R}$, and my received signal is $x[n]$, and my sample rate is $F_s$. I can get a signal $y[n]$ such that

$$y[n] = x[n]\cdot x[n-T_s F_s]$$

and I can probably decode everything from there, but it might be easier to sync if I used a decode signal, $d[n]$ such that

$$d[n] = d[n-1] + y[n]-y[n-T_s F_s]$$

Then $d[n]$ can be used to synchronize and to determine if the bit was a 1 or 0. I'm not some of the terms and variables I'm using aren't correct, and this might not be the best way to do it, but it's the only way I know how, would I would like to learn other methods of doing this, either in the time or frequency domains, to see if another method might fit one application better than another.


texts:
I need to calculate BER (bit error rate) for a system using spreading spectrum DS over ultra wideband (UWB). Modulation is bpsk. The spreading sequence lengths are not equal for all users. Receiver is matched filter. I know the equations of DS-UWB transmission signal and template signal which I should use at receiver to extract the user's data. I know that I should use convolution between the received signal (transmitted + AWGN noise) and template signal. My problem is, how to calculate the result of that convolution? is it integration or multiplication?
thanks in advance


texts: Orthogonal signal generator using integer arithmetic

texts:
I have a problem with implementing an orthogonal signal generator (OSG) algorithm on a microcontroller using integer arithmetic. I use this algorithm for a single-phase phase-locked loop (PLL) algorithm, for which I need an orthogonal component of a grid voltage.

The OSG algorithm is defined as follows:

$$\frac{d}{dt}v_x = \hat{\omega} \cdot \bigl((v_g-v_x)-v_y\bigr)$$
$$\frac{d}{dt}v_y = \hat{\omega} \cdot v_x$$

where $v_g$ is the measured grid voltage, $\hat{\omega}$ is the estimated grid frequency, and $v_x$ and $v_y$ are estimated components, with $v_x$ being equal to $v_g$ for ideal estimation. For this purpose, let us assume that the grid frequency is known.

The numerical integrator is implemented as follows:

$$y_k = \frac{T_s}{12} \bigr( 23u_{k-1} - 16u_{k-2} + 5u_{k-3} \bigl) + y_{k-1}$$

where $T_s=50~\mu\text{s}$ is the sample time.

Now, this algorithm works fine in floating point implementation, but is poor in integer arithmetic implementation. Here I give both implementations:

Floating point implementation

Code declaration.

float w = (2*PI)*50;
float Ts = 50e-6;
float i1u1, i1u2, i1u3, i1y1;
float i2u1, i2u2, i2u3, i2y1;

float NumInt3rd(float u1, float u2, float u3, float y1) {
    return (Ts/12)*(23*u1-16*u2+5*u3)+y1;
}

Main function.

float vg = floor(Input(0));

float vg_x = NumInt3rd(i1u1,i1u2,i1u3,i1y1);
float vg_y = NumInt3rd(i2u1,i2u2,i2u3,i2y1);

i1u3 = i1u2;
i1u2 = i1u1;
i1u1 = ((vg-vg_x)-vg_y)*w;
i1y1 = vg_x;

i2u3 = i2u2;
i2u2 = i2u1;
i2u1 = vg_x*w;
i2y1 = vg_y;

The Input(0) is a macro to get an input signal (sine wave with an amplitude of $2048$).

Integer arithmetic implementation

Code declaration.

int w = 643398L; // (2*PI)*50*2048
int i1u1, i1u2, i1u3, i1y1;
int i2u1, i2u2, i2u3, i2y1;

int NumInt3rd(int u1, int u2, int u3, int y1) {
    int iu = 23*u1-16*u2+5*u3;
    int iy = 240000L*y1;
    return (iu+iy)/240000L;
}

Main function.

int vg = (int) Input(0);

int vg_x = NumInt3rd(i1u1,i1u2,i1u3,i1y1);
int vg_y = NumInt3rd(i2u1,i2u2,i2u3,i2y1);

i1u3 = i1u2;
i1u2 = i1u1;
i1u1 = ((vg-vg_x)-vg_y)*w/2048;
i1y1 = vg_x;

i2u3 = i2u2;
i2u2 = i2u1;
i2u1 = vg_x*w/2048;
i2y1 = vg_y;

Note that I've checked for possible overflows, it never occurs. Also, the interesting thing is that the same algorithm works fine for $T_s=300~\mu\text{s}$.

I'm not that experienced with integer arithmetic implementations. Can you please give me an advice how to possibly fix this? Thanks!

Here is the estimation error for both implementations (on y-axis: percentage of the estimation error). The estimation error in case of integer artihmetic implementation is around $\pm 5\%$.


texts:
I am trying to calculate the change in pixels across two images.
Visually looking at them, there seem to be too many changes but calculating the distance between pixels gives me around about the same percent change values.

I am trying to calculate the distance here by adding the squared differences for each pixel and its square root which gives me the distance of value at rgb level for each pixel. i.e how much has the pixel value changed (sqrt((img1.r - img2.r)**2+(img1.g - img2.g)**2+(img1.b - img2.b)**2).

Then I am trying to check how many such pixels have changed by a threshold of 10.

I am not able to figure out, why does it show about the same difference in any two pair of images.

Is there a problem in the program or have I been thinking incorrectly

import pylab,skimage
import imageio,math,numpy as np

img1 = imageio.imread("029.jpg")
img2 = imageio.imread("028.jpg")
img3 = imageio.imread("030.jpg")
wh = img1.shape[0]*img1.shape[1]
squared_diff = ((img1 - img2)**2)
squared_diff1 = ((img3 - img2)**2)
squared_diff2 = ((img1 - img3)**2)
print float((np.sqrt(squared_diff[:,:,0] + squared_diff[:,:,1] + squared_diff[:,:,2])>10).sum())/wh*100
print float((np.sqrt(squared_diff1[:,:,0] + squared_diff1[:,:,1] + squared_diff1[:,:,2])>10).sum())/wh*100
print float((np.sqrt(squared_diff2[:,:,0] + squared_diff2[:,:,1] + squared_diff2[:,:,2])>10).sum())/wh*100

Output :

45.2615017361
46.1771918403
46.2858072917

For higher thresholds like 20 it gives 0 . np.max gives max distance as 15.x  everytime.

Is there any other way to find a feature like number of pixel changed with a threshold?

Any help is appreciated.


texts: MATLAB, recover a signal that has a cosinusoidal interference and that has been AM modulated

texts:
I have a signal which has been AM modulated but the channel also added a cosinusoidal interference A cos(2 * Pi * f0 * n * Tc + Phase), at the same f0 frequency of modulation. Tc is the sampling period.

I am given the distorted signal x(nTc) and I need to recover the original signal v(nTc)

I also have the filter coefficients of two filters:

A) a notch filter centered at 16KHz;

B) a low pass filter;

I started by calculating the fft of x(nTc) and it has a big peak at 7500 Hz, which is my f0
Now I need to eliminate the cosinusoidal interference and the notch filter seems to be the right way to do it but it is centered at 16KHz, so I have to center it at 7500Hz. To do so I used this Fourier transform property: exp(j * w0 * n)x(n) -> X(j(w-w0)). In my case I need to multiply the filter coefficients of the notch filter by exp(j * 2 * Pi * 8500 * n * Tc)
Everything should be good, so now I can filter the x(nTc) signal using input_without_interference = filtfilt(notch_coeff,1,x(nTc))

At this point I have to AM demodulate the signal. Demodulation is accomplished by multiplying the signal by the same carrier that modulated it, in this case cos(2 * Pi * f0 * n * Tc)
To do so I used this code:

for k=0:length(input_without_interference))-1
    dem_signal(k+1) = cos(2 * Pi * 7500 * k * Tc) * input_without_interference(k+1);
end

Now If I plot the signal dem_signal it still shows some traces of modulation so something is certainly wrong.

Assuming everything went good, I now have many replicas of the original signal v(nTc) and one of them is centered in the origin, so to take it I just need to use the low pass filter. v(nTc) = filtfilt(lowpass,1,dem_signal)

Result: the output is still distorted and it has an imaginary component! Where is my
 mistake?


texts: Converting speech audio to telephone audio

texts: I am trying to train a machine learning algorithm on telephony speech audio. However, there isn't really enough data for this anywhere that I can find. My solution is to just use speech audio from Youtube videos. How can I convert audio from a Youtube video to audio that simulates telephony audio?


texts: How to choose between rising and falling edge (e.g. when using a pulse generator)?

texts:
In StateFlow, one often uses the command e.g. on clk: count++; I wonder how the clock is counted exactly. Is it the rising edge, or the falling or does it count both?
When I enter the properties of the clock I can't see any setting about that and when simply using "on clk" I can't see it neither.


texts: Samples of period to frequency, sampling theorem

texts:
In the sample rate of 96kHz there are 96000 samples in a second.

So a periodic waveform with period 96000 samples has a frequency of 1Hz.

A signal with period of 9 samples is ~ 10666 Hz.

Of 5 samples 19200Hz

4 -> 24000

3 ->32000

2 -> 48000 (Nyquist)

So the question is for high frequencies are there really so big gaps in sampling. It can't really represent a frequency of 28.500?

Edited:

I made some tests creating waveforms in Wavelab. When i try to create a waveform with frequency ( sample rate / 2 ) I get nothing. 0 samples.

But when i create any other frequency ( < Nyquist ) like 47999 in 96kHz sample rate a get visually a kind of amplitude modulated waveform that in the FFT window  it shows the correct frequency.

So as i understand the frequency is represented as amplitude modulation of successive samples (how this works i have to find out). Funny thing is that for the waveform to actually resemble the sought form the frequency must be ~ 1/10 of the sample rate.


texts: Mathematical proof and intuitive reasoning for: $\lim \limits_{\epsilon \to 0} \frac{1}{\pi t}\sin\frac{\pi t}{\epsilon}=\delta(t)$

texts:
Please provide the mathematical proof and intuitive reasoning for:

$$\lim \limits_{\epsilon \to 0} \frac{1}{\pi t}\sin\frac{\pi t}{\epsilon}=\delta(t)$$
Where $\delta(t)=$Unit impulse function

Intuitive reasoning (what I have gathered):

  If we compare $\sin\frac{\pi t}{\epsilon}$ with $\sin(\omega t)$, then we
  can see that with decreasing value of $\epsilon$ the waveform $\sin\frac{\pi
 t}{\epsilon}$ approaches a sine wave of very high frequency.

  Along with that, $\frac{1}{\pi t}$ as a multiplicative factor implies
  that at $t=0$, the function would be infinity and as the value of $t$
  increases, the amplitude of  $\sin\frac{\pi t}{\epsilon}$ goes on
  decreasing with the function tending to zero as t tends to infinity.

Is my reasoning correct?

Please help me with the mathematical proof as well.

Thank you :).


texts: Why 8-PSK modulation has higher Eb/N0 than QPSK or BPSK?

texts: I run a transmitter/receiver simulation using QPSK, BPSK, and 8-PSK modulations. As far as I have noticed, the best SNR per bit results are obtained when simulating with 8-PSK modulation (and also at longer distances). For example - at a distance of 10 m (transmitter to receiver) Eb/N0 for QPSK was 22.9 and for 8-PSK it was 26. Does it depend on transmitting more bits symbols per cycle (8-PSK transmits 3 while QPSK/BPSK only 2)?


texts: What is the ROC for this discrete signal:

texts:
$$ x(k)=4[u(k-2)-u(k)*δ(k-3)]$$

I found that the $\mathcal{Z}$ transform of the signal is $X(z)=4/(z^2)$.

What would the ROC be?


texts:
The super resolution model is as such: Y_k = DHFX + N , where Y_k is one of the N measured low resolution image, D is the Decimation Matrix, H is the blur matrix , F is the warp matrix, N is the noise and X is the high resolution image. Lets say the X is [L *L] or we can say [L^2 * 1], F is [L^2*L^2], and H is [L^2*L^2] and D is [M_k^2*L^2] where M_k is the size of Y.
In the super resolution problem we need to find X but for that we need to find the F(warp matrix). The shifting is done to get the sub pixel information as far as I understand. I am also able to find the [3*3 ] warp matrix which transform the coordinate space but that is not required. So if anyone can help me to find the F matrix. Thank you


texts: How to find the phase spectrum of a rectangular pulse? (Fourier Transform)

texts:
How to find the phase spectrum of a rectangular pulse?

The Fourier transform of a rectangular pulse

$$ x(t) =
\begin{cases}
1,  & \text{for $|t| \le \tau /2$ } \\
0, & \text{otherwise}
\end{cases}
$$

is given by:

$$F[x(t)]=\tau [\frac{\sin \omega(\tau /2)}{\omega (\tau /2)}]$$

In general, the Fourier Transform, $X(\omega)$ is id a complex valued function of $\omega$. Therefore, $X(\omega)$ can be written as:

$$X(\omega)=X_R(\omega)+jX_I(\omega)$$

The magnitude of $X(\omega)$ is given by

$$\vert X(\omega) \vert= \sqrt {(X_R(\omega))^2+(X_I(\omega))^2}$$

The phase of $X(\omega)$ is given by

$$\angle{X(\omega)}=\tan^{-1}\frac{X_I(\omega)}{X_R(\omega}$$

  Question:

  1)How can we find the phase spectrum of a rectangular pulse as there doesn't appear to be any imaginary part in $X(\omega)$?

  2)Why is the phase spectrum changing the way it is when we are time
  shifting the rectangular pulse from the origin?

  3) Any other insight regarding the magnitude and phase spectrum is
  welcome.


texts: How to detect peaks ignoring soft declivity?

texts:
I have to detect repetitions in a signal, and I am managing to do it by finding peaks and valleys, and cutting the desired pieces. The following figure shows the pattern wanted, shaded areas are the portions of interest:

The problem is: some of the signals that I'm analyzing, instead of starting with a constant zero, show a very smooth decline, as shown in the following figure:

As a consequence, a peak is detected as one of the first values, and the first repetition is detected wrongly. I'm looking for suggestions to avoid it. I have already tried (not sure if correctly):

Using derivatives (diff) to detect inclination.
Using a robust regression (robustfit) and finding inclination by the coefficients - too slow.

I'm using a MATLAB code created by someone else, modified by me (original can be found here):

function repeticoes = peakdet(v, delta, x)
%PEAKDET Detect peaks in a vector
%        [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
%        maxima and minima ("peaks") in the vector V.
%        MAXTAB and MINTAB consists of two columns. Column 1
%        contains indices in V, and column 2 the found values.
%
%        With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
%        in MAXTAB and MINTAB are replaced with the corresponding
%        X-values.
%
%        A point is considered a maximum peak if it has the maximal
%        value, and was preceded (to the left) by a value lower by
%        DELTA.
% Eli Billauer, 3.4.05 (Explicitly not copyrighted).
% This function is released to the public domain; Any use is allowed.

maxtab = [];
mintab = [];

v = v(:); % Just in case this wasn't a proper vector

if nargin < 3
  x = (1:length(v))';
else
  x = x(:);
  if length(v)~= length(x)
    error('Input vectors v and x must have same length');
  end
end

if (length(delta(:)))>1
  error('Input argument DELTA must be a scalar');
end

if delta <= 0
  error('Input argument DELTA must be positive');
end

mn = Inf; mx = -Inf;
mnpos = NaN; mxpos = NaN;

lookformax = 1;

for i=1:length(v)-1
  this = v(i);

  if isempty(maxtab)
    if this >= mx, mx = this; mxpos = x(i); end
  else
    if this >= mx, mx = this; mxpos = x(i); end
    if this < mn, mn = this; mnpos = x(i); end
  end

  if lookformax
    if this < mx-delta
      maxtab = [maxtab ; mxpos mx];
      mn = this; mnpos = x(i);
      lookformax = 0;
    end
  else
    if this > mn+delta
      mintab = [mintab ; mnpos mn];
      mx = this; mxpos = x(i);
      lookformax = 1;
    end
  end
end
    if size(maxtab, 1) == 0
        inicios = [];
        finais = [];

    elseif size(maxtab,1) >= size(mintab,1)
        inicios = maxtab(1:size(mintab,1),1);
        finais = mintab(:,1);
    elseif size(maxtab,1) < size(mintab,1)
         inicios = maxtab(:,1);
        finais = mintab(1:size(maxtab, 1),1);
    else
        disp('Erro ao reconhecer repeticoes!')
    end

    repeticoes = [inicios finais];

end

Any suggestions are welcome!


texts:
I have read that we can use audio watermarks to "identify the original produces of the material or the consumers of the media."

How consumers of the media is measured by audio watermarks?

Lets say in TV we play an video which has audio watermark in it and 100's of people are watching it. How do we capture the audio watermark that got played in viewers TV?


texts: use chirp z transform for spectral resolution

texts:
let us suppose that we some signal $x(t)$,which consist of sinusoidal   components and white noise,i would like to know how to use    chirp z transform for improving spectral resolution?i found   this article about  chirp z transform

http://prod.sandia.gov/techlib/access-control.cgi/2005/057084.pdf

also  similar  article on official matlab side

http://www.mathworks.com/help/signal/ref/czt.html

but  as i know always in spectral estimation methods,some  knowledge we should know how process is going,or  some other things,let us suppose that we know sampling frequency lat say $Fs=100$,and   also we know that it consist as i mentioned by periodic components(sin or cosine),but i dont  know number of components,neither frequency or neither phases(+neither amplitudes),how can i apply chirp z transform to this signal?we have just some signal samples  with size $N=294$ ,please help me how to  choose correct components so that  get  effectively spectral structure and if let say one frequency is $f_1=13$ and second $f_2=13.5$,this transform be able to separate this  frequencies from each other,thanks


texts: Combining spectrograms with different windows to getting arbitrary time and frequency resolution

texts:
Spectrograms convert an amplitude-time signal into a frequency-time signal using Short Time Fourier Transform (STFT). The window size L used in STFT is a design choice. As L increases, you get higher frequency resolution and lower time resolution and vice versa.

If you performed STFT with a large L to get spectrogram A with high frequency resolution then performed STFT with small L to get spectrogram B with high temporal resolution, can you combine spectrogram A and B somehow to get high temporal and frequency resolution?

My gut feeling is no, as changing L is simply "changing how you interpret the amplitude-time data."


texts: How can I identify Arabic letters in an image using Java?

texts:
I want to identify Egyptian car number plates that look like this:

How can I do this in Java ?

I plan to do this in Android, as the user will catch the number plate with a camera and the results will be submitted to server.


texts: Is time varying channels same as circulary symmetric complex gaussian?

texts: From the literature i understood that if the channel is non time varying then we can consider it as circularly symmetric complex gaussian distributed. So can we consider the same if the channel is time varying ?


texts: How to detect malformed-ellipse-alike objects?

texts:
I need to process a stack of ultrasonic images, where there is mostly a lot of noise to see, but I need to extract specific objects that look for example like this:

Here is an example from Wikipedia of  pretty nice formed object

But I'm not sure how to get the whole region.

What I did:

I already tried with mean & median filters to achieve a binary image in the end, which then is processed with the union-find-method to create regions. Unfortunately this leads to vague results. And it seems like the difference of the edge pixels aren't high enough to get a clean sobel-image (for trying to define the margins).

What I believe could work

On a different project I was working with the hough-algorithm to detect a (very clean shaped) ball in a picture. I heard this method can be applied for lines and thus for curves too. But would I be able to detect a region like this ? And if so, how would it be? is it multiple hough-detection for just one area? Since I need to process in real time I'm afraid this would exceed the processing frame (or wouldn't it?)

My question

Has anyone an idea how to easily detect such objects ? Would it be possible to apply hough under real-time conditions, on an ordinary (not an high end) computer too?

Thanks for any advises!
Kind Regards


texts:
Today, in my lecture on intersymbol interference (ISI), there was a discussion on postcursor ISI, and precursor ISI. Postcursor ISI is caused by the past bits, whereas, precursor ISI is caused by the future bits on the present bit. I could not understand how precursor ISI can occur. I was told that it is because of the delay of the channel. What does it mean?
Isn't there a causality violation? What am I missing?

Can someone please explain? Thanks.


texts:
Can you tell me whether the following signal is a energy or power signal?
$$
x(n) = e^{j (n \pi/2 + \pi/8) }
$$
I've solved it and found it as it was neither power nor energy signal was that right?


texts: The effect of upsampling on DFT coefficients

texts:
I am learning DSP by myself and I encountered a problem that bewilders me.
If I have a sequence of length N, and I upsample it by a factor of 3. How would the DFT change or related?

For example:

x = [1 2 3 4 5 6 7];

y = [1 0 0 2 0 0 3 0 0 4 0 0 5 0 0 6 0 0 7];

What is the difference between $X[k]$ and $Y[k]$?


texts: Threshold for peak detection of noisy signal in frequency domain

texts:
I have a noisy signal (Gaussian Noise) with a known SNR and known Noise variance.
This signal contains a number of Frequency components.
My ultimate objective is to detect the components of the signal.
What is the Threshold formula should I use to detect the peaks of the signal in frequency domain given a certain Pfa (probability of false alarm)?

The following plot is an example of my objective. I want a Threshold that separates my useful components from the noise components (Let's say between 0.6 and 1)


texts: When does the Amplitude Spectrum equal the Magnitude Spectrum?

texts:
I had the following question on edX:

I'm failing to understand why the second signal has $M(\omega)=A(\omega)$. First I find the DTFT of the signal:

$$\mathcal{F}\{\delta [n]+\delta [n-1]\}\ =\ 1\ +\ e^{-j\omega}$$

I then deduce the $A(\omega)e^{\ j\phi_A(\omega)}$ representation:

$$1\ +\ e^{-j\omega}\ =\ \big[\cos (\omega)+1\big ]\ -\ j\cdot\sin (\omega)\ =\ 2\cos \bigg (\dfrac{\omega}{2}\bigg )e^{-j\frac{\omega}{2}}$$

The magnitude spectrum is:

$$M(\omega)\ =\ 2\cdot \bigg |\cos \bigg (\dfrac{\omega}{2}\bigg )\bigg |$$

Therefore $M(\omega)\neq A(\omega)$. As a result, $\phi_M(\omega)$ will equal $\phi_A(\omega)$ except contain discontinuities of size $\pm\pi$ at values of $\omega$ where $A(\omega)$ has a zero crossing.

The only scenario I can envisage where $M(\omega)=A(\omega)$ is if $A(\omega)\geq 0$ however as I got the question wrong there must be something wrong with my logic, could somebody clarify?


texts: How do I stabilize my oscillating system?

texts:
I have two scalars $x$ and $y$ that vary with time $t$ such that

\begin{align}
\frac{dx}{dt} &= ay + b + dx\\
\frac{dy}{dt} &= cx
\end{align}
 but $a, b, c$ are unknown.

If $d$ is too small, the system will be underdamped (and oscillate badly), but if $d$ is too large, the system is overdamped, and doesn't converge fast enough.

How do I set $d$ based on statistics taken about historic $x, y$ so that the system is critically damped and moves as quickly as possible to the point where $x=0, y=\bar y = -\frac{b}{a}$?


texts: Adaptive Thresholding for liver segmentation using Matlab

texts:
I need to segment liver from an abdominal CT image using Adaptive Thresholding. But I get the entire foreground separated from the background alone. I need only the liver part of the foreground separated.
Check out the pdf file in http://www.ijcaonline.org/casct/number1/SPE34T.pdf
I need an output similar to that shown in Figure 6.

I attach my coding here. Kindly help me out.

%testadaptivethresh.m
clear;close all;
im1=imread('nfliver2.jpg');
bwim1=adaptivethreshold(im1,11,0.03,0);
figure,imshow(im1);
figure,imshow(bwim1);
imwrite(bwim1,'at2.jpg');

function bw=adaptivethreshold(IM,ws,C,tm)
%ADAPTIVETHRESHOLD An adaptive thresholding algorithm that seperates the
%foreground from the background with nonuniform illumination.
%  bw=adaptivethreshold(IM,ws,C) outputs a binary image bw with the local
%   threshold mean-C or median-C to the image IM.
%  ws is the local window size.
%  tm is 0 or 1, a switch between mean and median. tm=0 mean(default); tm=1 median.
%
%  Contributed by Guanglei Xiong (xgl99@mails.tsinghua.edu.cn)
%  at Tsinghua University, Beijing, China.
%
%  For more information, please see
%  http://homepages.inf.ed.ac.uk/rbf/HIPR2/adpthrsh.htm

if (nargin<3)
    error('You must provide the image IM, the window size ws, and C.');
elseif (nargin==3)
    tm=0;
elseif (tm~=0 && tm~=1)
    error('tm must be 0 or 1.');
end

IM=mat2gray(IM);

if tm==0
    mIM=imfilter(IM,fspecial('average',ws),'replicate');
else
    mIM=medfilt2(IM,[ws ws]);
end
sIM=mIM-IM-C;
bw=im2bw(sIM,0);
bw=imcomplement(bw);

My modified code for testadaptivethresh.m

clear;
im=imread('nfliver7.gif');
figure,imshow(im)
bwim1=adaptivethreshold(im,300,-0.15,0);
bw=bwareaopen(bwim1,3000);
se=strel('diamond',4);
er=imerode(bw,se);
bw1=bwareaopen(er,3000);
er1=imerode(bw1,se);
bw2=bwareaopen(er1,1000);
fi=imfill(bw2,'holes');
figure,imshow(fi)

op=uint8(fi);
seg=im.*op;
figure,imshow(seg)
imwrite(seg,'sliver7.jpg');


texts: What are the advantages and disadvantages of the Ridgelet transform?

texts:
I used the continuous wavelet transform (CWT) to analyze the behavior of a 1-D smooth signal (extracting its modal property, damping and frequency components) in my research. Now someone asked:

  Why the presented CWT algorithm is better than Ridgelet transform?

I thought it is better to ask in this forum about disadvantages of Ridgelets.
For example:

it is said that Wavelet can handle point discontinuities whereas
Ridgelet is better to deal with line discontinuities. So, do we have
any line discontinuities in a 1-D signal?
what is computational complexity of Ridgelet in term of big $O$ i.e.
(like $O(N)$), for both continuous (in scale) and discrete versions?


texts: Image Restoration and Standard Forms of Second Order Cone Programming (SOCP)

texts:
I'm studying the application of SOCP methods in Image restoration  And I want to understand the difference between the two formulas of SOCP and how they are related.
Standard form (1) :

min $f^{t}x $

subject to $\left \| A_{i}x + b_{i} \right \|_{{2}} \leq c^{{i}^{t}}x + d_{i}$

but also in some papers i found this standard form (2) :

min $c{_{1}}^{t}x_{1} + c{_{2}}^{t}x_{2} +. . . +c{_{r}}^{t}x_{r} $

$A_{1}x_{1}+ A_{2}x_{2} +...+A_{r}x_{r}=b$

What are the major differences between the two?
when and how should each one be applied?


texts: Methods for calculating signal quality metric

texts:
i'm sure this is a fairly large topic area, but I'm not quite sure where to start. I have a signal in which I can change several parameters, such as the frequency, gain, etc. The modulation is a differential BPSK, and my current metric for determining the signal quality is what percentage of my auto correlation (correlating the current bit signal with the previous bit) is above 0. If it's all above or below zero, a 100% signal quality. if it's half and half, then it's 0%, and then it's linear from there. Is there a better way of doing this? I'm wanting to be able to change my signal parameters in order to increase my signal quality, where a higher signal quality should be a higher reliability of receiving a message correctly.

I'm also looking at doing demodulation in the frequency domain, so then my calculations would probably need to change. Any input and references would be appreciated. Thanks!


texts: Discrete Time Fourier Transform to Frequency Domain

texts: I failed to do this question on the exam and finding it very difficult, I would be glad if you can help me solve it. How shall I start?


texts: If you shine a bright enough white light on a red object, will the object turn white or stay red? Please assume the object is not capable of specular reflection.


texts: Discrete State Space Model - Why Are We Calculating $ x \left[ k + 1 \right] $ Instead of $ \dot{\boldsymbol{x}} \left( t \right) $?

texts:
A continuous state space model is defined as follows.
$$
\dot{\boldsymbol{x}}(t) = A \boldsymbol{x}(t)+ B \boldsymbol{u}(t) \\
\boldsymbol{y}(t)= C \boldsymbol{x}(t)+ D \boldsymbol{u}(t)
$$

If we want to calculate the state $\boldsymbol{x}$ for a specific time $t$, we need to use the state transition matrix $\Phi$. However, when looking at a discrete state space model, things look different.
$$
\boldsymbol{x}[k+1] = {A}_{d} \boldsymbol{x}[k] + {B}_{d} \boldsymbol{u}[k] \\
\boldsymbol{y}[k] = C \boldsymbol{x}[k] + D \boldsymbol{u}[k]
$$

With the discrete state space model, there seem to be no 1st order differential equations incorporated anymore. Instead, we can directly calculate the next sample, without the need for a state transition matrix. How come this is the case?


texts: Filterpy Kalman Filter batch processing with multiple measurement sources

texts:
In pythons module for kalman-filtering, filterpy, there is a function batch_filter() to batch filter a list of measurements that then can be used for RTS-smoothing. See the documentation here.

I want to smooth my measurements with a RTS-smoother like:

(mu, cov, _, _) = kf.batch_filter(list(np.array(centroids)))
(x, P, K, Pp) = kf.rts_smoother(mu, cov)

The problem is now that I have two measurements from two different sources with different measurement noises. The function batch_filter() can only process one source of measurements. When I calculate and save x and P for every time step with the sequence "predict, update with R1, update with R2", is this the same thing batch_filter() would do?


texts:
I am calculating the STFT using scipy, and then reducing it in size, flattening and then normalising it. Code is below. It's incredibly slow, given I'm trying to process 64,000 audio streams of about 10 seconds each. Does anyone have any suggestions on how I might speed this up?
f, t, Zxx = signal.stft(data, fs=sample_rate, nperseg=256)
def normalise(z):
    return np.log(np.abs(z + 1e-12)) + 1e-7

# Shrink size of STFT
avg_window = 256
a = np.zeros((Zxx.shape[0], int(Zxx.shape[1]/avg_window)))
for j, row in enumerate(Zxx):
    for i in range(int(Zxx.shape[1] / avg_window)): # note: truncated at end
        a[j][i] = np.mean(np.abs(row[i*avg_window:(i+1)*avg_window]))

flattened = np.array(Zxx).flatten() # 129 x 54
normalised = np.array([normalise(z) for z in flattened])


texts: Object detection versus object recognition

texts:
I have a slight confusion differentiating between object recognition and object detection. Some people say object detection is a sub-topic of object recognition? Can someone clarify the the difference between these two topics?

To the best of my knowledge:

Object Recognition is responding to the question "What is the object in the image" Whereas, Object detection is answering the question "Where is that object"?

Hope someone can illustrate the difference by also generously providing an example for each.


texts: Combining CDMA with differential Manchester-like encoding?

texts:
I wonder to which extend it is possible to combine a self-clocking code (like Differential Manchester encoding) with synchronous CDMA, and how it would be done?

Suppose I want to receive some bits from devices where each of them has sufficiently synchronized clocks but the receiver of the bits does not (and hence suffers from clock drift, delays etc).

Right now (without CDMA) I do not even need to care about the clock ... I just detect the transitions at the receiver and can infer the bits (because there is a guaranteed transition every bit). I wonder if something similar can be done while still using CDMA?

Take for example 3 IoT devices which send data in parallel via 4-order Walsh-Haddamard code.

The code matrix is:

$$
\begin{bmatrix}
  1 & 1 & 0 & 0 \\
  1 & 0 & 0 & 1 \\
  1 & 0 & 1 & 0 \\
  1 & 1 & 1 & 1 \\
\end{bmatrix}
$$

and the three devices send bits $001011010101$, $001010110011$ and $001010001111$. The waveform as seen at the receiver looks like:

(as per assumption above they sum up nicely in the channel and preserve orthogonality of the data bits but the receiver does not know the clock). We can see that there is not necessarily a transition (or at least change of sample value) every bit/symbol. I tried replacing each 1 with $[1, 0]$ and each 0 with $[0, 1]$ in the code matrix but this doesn't change anything significantly.

What would be the best way to implement this?

(I hope my question is clear enough. If not, please leave a comment and I will clarify)


texts: Complex Distance Transform containing distance and angle info?

texts:
Is it possible to obtain a distance transform with the distance (float) and angle (in radians) from the nearest pixel?
I would be using python and the library I am currently relying on is opencv-python, although I am open to whatever works.


texts: If we define samples manually how could there be a Sampling Frequency?

texts:
I just got a thought when generating matrices , plotting it and taking its FFT,

Lets say, I have randomly input a matrix of length 256,

[1 7 9 5 7 1 -2 0..................8]

According to the FFT, the resulting 256 FFT coefficients correspond to the following frequencies:

0, $f_s/256$, $2 f_s/256$, $3 f_s/256$, $\ldots$, $255 f_s/256$

but since the samples i taken has no sampling frequency defined how could there be some sampling frequency ? Like in an image ,

Since the image is just a 2 dimensional Matrix , if i randomly generate an image in a matlab with random pixel values (sample points) and plot its FFT , the matlab will show the frequencies of that 2-D matrix but since there is no Sampling Frequency is defined how matlab calculates the frequency for each of the coefficient ?


texts: Conjugate symmetry of the DFT of real-valued sequences

texts:
I have read about Fourier transformation that real signals are "mirrored" in the real and negative halves of the Fourier transform because of the nature of the Fourier transform.  For example, if we have the real signal $x$ whose length is $N$, its $X = FFT(x)$ should be conjugate symmetric. That means the $X^*[m] = X[N-m]$.
That is ok, and I have read it clearly HERE
My question: assuming you have a real signal $x = [-1, 1, -1, 1];$ its $X = FFT(x) = [0,     0,    -4 ,    0]$
So, in that case, we cannot say that $X^*[m] = X[N-m]$.  Why ?
Is what I said above is right ? in which cases is right?


texts: Comparing signals recorded using diffrent methods

texts:
Suppose we have two signals( focus to audio) recorded with different methods,
or more simply, suppose we have two audio signals recorded with different microphone, and  we want to confirm that these signals are of the same audio.

What are the available alogorithms or libraries that make this possible?


texts: algorithm for feature extraction from audio signals

texts: I am trying to build a model for speaker identification, and I understand that the first step is to extract the features from the audio signals that are in my database. I have done quite a bit of research and cant find how to do this extraction and to which features. How do I select which features to extract to? what would be the algorithm? Can someone please point me to either an algorithm that I could implement or a sample matlab code?


texts:
I have been given the task of doing audio compression using FFT.
The assignment is below.

Below is my MATLAB code:

% part 1
[y, f] = audioread('file.wav');
y = y(1:1024*161);
% p = audioplayer(y,f);
% play(p)
Ya = fft(y);
Y = fftshift(abs(Ya));
N = length(Y);
w = -f/2+f/N:f/N:f/2;
figure(1)
plot(w,Y)
% part 2
Y1 = Ya((N/4):(3*N/4));
% part 3
y = ifft(Y1);
player = audioplayer(y,f);
play(player)

Now here is my problem. I was able to do part 1 and part 2 perfectly with no problem, but for part 3 the problem begins. I get no sound at all!!. Now in the prompt it asks to do FFTSHIFT before the IFFT, to me that does not matter, either how I get now sound. can somebody help me on how to fix this please?


texts: Are there any State Space Control C libraries?

texts:
I'm searching for a C library which implements methods needed for state space control (e.g. state feedback controller and state observer equations).
Googling various terms only brought up a non-C (C# ?) library:
https://github.com/tomstewart89/StateSpaceControl
Therefore I was hoping to find existing projects which use C state space control algorithms or ready to use libraries with your help.
Thanks in advance!


texts:
I'm studying what I believe to be a system that produces a (noisy) quantal response (i.e. there should be a minimum increment of response). To do this, I've plotted a histogram of some ~6000 amplitude values from the system into small bins, and performed an FFT on the histogram, expecting that a common quantal value should appear as the lowest frequency "spike". However, my spectral response:
 looks odd to me. What I see are what look like harmonics.

Why does this occur in the quantal release problem? Is this lowest frequency actually my quantal response, or is it a typical electrical response/artefact from the FFT of a signal with so many discontinuities? I thought that odd-numbered harmonics are the most powerful, but here they seem to be linearly declining. Why am I seeing other smaller spikes appear 1/3rd of the way through each of the "main" harmonic bands?


texts: What plotting / visualization tools are useful in generating good plots?

texts:
It has been said a picture is worth a thousand words!
There have been a large number of softwares to write dsp codes and analyze signals and algorithms. It will be good to share on data plotting tools, that is - given a data set (1-d, 2-d, 3-d, ..higher dimensions (creative ways to visualize this)) what are the different tools available and advantages of each for plotting such data sets.
Advantages can be with respect to:
(a) quick/easy plotting
(b) formats of fugures (.eps, .svg) file which can be generated
(c) latex font embedding easiness
(d) post script figures
Tools like: matlab, octave, gnuplot, tikz, pgf etc are widely used by many of us. It would be useful to share your experiences (advantages/disadvantage) and best plot you have drwan using the tool you prefer!


texts: Ratio of the power in the sidebands to the total power in conventional AM

texts:
Message signal and carrier signal are given as

$m=\cos(2 \pi  \cdot f_{m} \cdot t) $

$c=A_{c}\cdot \cos(2 \pi  \cdot f_{c} \cdot t) $

I was asked to find modulation index $a$ so that ratio of the power in the sidebands to the total power in the transmitted signal will be $\frac{1}{3}$

My attempt

We know that, in conventional AM, total power of transmitted signal $P_{u}= \frac{Ac^2}{2} + \frac{Ac^2}{2}\cdot a^2 P_{m} $  where $P_{m} $ is message signal's power.

I think that power in sidebands is related to $\frac{Ac^2}{2}\cdot a^2 P_{m}$ term

$\frac{Ac^2}{2}\cdot a^2 P_{m} = \frac{1}{3}P_{u}$

Submitting $P_{u}$

$\frac{Ac^2}{2}\cdot a^2 P_{m} = \frac{1}{3}(\frac{Ac^2}{2}(1 +  a^2 P_{m}))$

$ a^2 P_{m} = \frac{1}{3}(1 +  a^2 P_{m})$

$ a^2 =  \frac{1}{2P_{m}}$

It seems false because modulation index might have been a numerical value.


texts: Evaluating different inputs into a given system

texts:
We need to find which signal wont be the same after going through this sytem.

My approach is this one:

I took the filters frequency and did 2pi * 1/4 and gave me this "limit" 1.5707

(2pi because the periodicity and 1/4 because the filter "bandwith"?)

Then i proceed with the inputs evaluating them this way:

a) 4 * pi/7 = 1.7951
b) 4 * pi/9 = 1.3962
c) 4 * pi/8 = 1.5707
d) idk

The result "a", which is out of the "limit" is the one that does not follow the rule y[n]=x[n] and thus, the correct answer. (It is the actual correct answer)

Did i have luck there with the math or the approach is correct?


texts: Difference between $\tt rand()$ and $\tt randn()$ in MATLAB

texts:
rand() is for "Uniformly distributed pseudorandom numbers"

I generate two rand arrays, and then use the cross-correlation (xcorr2) and the normalized cross correlation (normxcorr2) as follows:
a = rand(512,512);
b = rand(512,512);
c = xcorr2(a,b);
figure;surf(c);shading flat;axis tight;
d = normxcorr2(a,b);
figure;surf(d);shading flat;axis tight;

The results are as follows:

in which there is strong correlation in the xcorr2().
randn() is for "Normally distributed pseudorandom numbers"

The same as the above:
a = randn(512,512);
b = randn(512,512);
c = xcorr2(a,b);
figure;surf(c);shading flat;axis tight;
d = normxcorr2(a,b);
figure;surf(d);shading flat;axis tight;

and the results are:

Then, my question is why there is apparent difference of xcorr2() for the two random data?


texts: How to make the harmonic complex tones sound comfortable?

texts:
I'm trying to create a complex tone, which is consisted by several sine waves with a 220Hz fundamental frequency and 6 multiplied harmonics (at 440Hz, 660Hz, 880Hz, 1100Hz, 1320Hz, 1540Hz).
My code below generated a tone. However, it sounds very uncomfortable, and its waveform and spectrogram picture also look much intensive.
fs = 44100;      % [Hz]
T = 0.5;          % [s]
x = 0:1/fs:T;   % [s]

f = 220;        % [Hz]
y = 1*sin(2*pi*f*x);
f1 = 0.2*sin(2*pi*440*x);
f2 = 0.2*sin(2*pi*660*x);
f3 = 0.2*sin(2*pi*880*x);
f4 = 0.2*sin(2*pi*1100*x);
f5 = 0.2*sin(2*pi*1320*x);
y1 = y + f1 + f2 + f3 + f4 + f5
fname = 'sin_220.wav';
audiowrite(fname,y1,fs);

I found that one research group have shared similar tone materials, in which, the sine waves were synthesized regularly and sound pleasant.

My question is how to make adjustments to the code so that it looks like this.


texts:
Analog AGC is mainly used to increase the signal level to match the ADC sensitivity? Is Analog AGC also used to attenuate the signal?.
Digital AGC is the first block used after ADC. Is the main purpose of digital AGC is to provide constant level to the signal processing receiver chain( RRC, CIC,  demod....)
So analog AGC output is range of signal levels but digital AGC always gives out a fixed level?


texts: Issue with demonstrating amplitude modulation of sinusoid using MatLab

texts:
I am trying to demonstrate Conventional and DSB-SC modulation of a signal $ \sin(2\pi1000t) $ using MatLab. My carrier signal is $ \sin(2\pi f_ct)$ where $f_c = 100\text{ MHz}$. Now according to theory, after modulation, the modulated signal should have a spectrum same as that of message signal shifted to the frequency of carrier on both sides. But I am getting the plot as that of the carrier signal in the frequency domain.
This is the code that I am using for DSB-SC modulation
% CODE FOR PLOTTING MESSAGE SIGNAL %%%%%%%%%%%%%%%%%%%%%%%%%%%%

Fs = 20*100*10^6/3;
Fm = 1000; %Hz
t = 0:1/Fs:0.005-1/Fs;

Am = 1;

m = Am*sin(2*pi*Fm*t);

temp = size(m);
N = temp(2);

f = -Fs/2:Fs/(N-1):Fs/2;

M = fftshift(fft(m));
figure(1)
plot(f,abs(M))
grid('on')
title('Message signal')
xlabel('Frequency in Hertz')
ylabel('|M(F)|')

% CODE FOR PLOTTING CARRIER SIGNAL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Fc=100*10^6;%100 MHz

Fs = 20*Fc/3;
Ac = 1;%amplitude
c = Ac*sin(2*pi*Fc*t);

temp = size(c);
N = temp(2);

f = -Fs/2:Fs/(N-1):Fs/2;

C = fftshift(fft(c));
figure(2)
plot(f,abs(C))
grid('on')
title('Carrier signal')
xlabel('Frequency in Hertz')
ylabel('|C(F)|')

% CODE FOR PLOTTING DSB-SC MODUALTED SIGNAL %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s = Ac*m.*sin(2*pi*Fc*t);

temp = size(s);
N = temp(2);

f = -Fs/2:Fs/(N-1):Fs/2;

S = fftshift(fft(s));
figure(3)
plot(f,abs(S))
grid('on')
title('DSB-SB Modulated signal')
xlabel('Frequency in Hertz')
ylabel('|S(F)|')

And this is the spectrum of $S(f)$ in code. As you can see. there are only two peaks, one at $+f_c$ and one at $-f_c$(not shown in figure). But in reality, it should have been 2 peaks at $f_c+f_m , f_c-f_m$ and two peaks at $-f_c-f_m,-f_c+f_m$. So what is going wrong here? How can I fix this?


texts: Image Contrast Enhancement based Sub-histogram Equalization Technique without Over-equalization Noise

texts:
I'm a computer science student and I have an assignment requiring me to implement the algorithm proposed by this paper but I'm hitting a roadblock trying to understand the implementation of one of the methods proposed in the paper. I'm not sure if I understand Section IIIB well.
This is what I understand so far.

In order to separate the histograms into subregions, we must first obtain the histogram of the intensity levels. Once achieved, we calculate the density of the histogram (equation 7).
Furthermore we have to sort the density function values in descending order. How is this implemented? I assume its only in the following paragraph in the paper that goes on to describe this but I'm having a hard time understanding what the authors are trying to convey.

Once we find the groups of high density areas, we extend a group area surrounding these values. The paper talks about momentary changes. I'm not sure what this means, anyone have an idea?

The page is located at
http://www.waset.org/journals/ijeee/v3/v3-6-53.pdf


texts: Proof of transfer function factorization $\frac{b_0}{a_0} \frac{\prod_{k=1}^M (1-c_kz^{-1})}{\prod_{k=1}^N(1-d_kz^{-1})}$

texts:
This is from Oppenheim's Discrete-Time Signal Processing, but the book doesn't seem to describe how the factorization is done.

The transfer function:

$$H(z)=\frac{Y(z)}{X(z)} = \frac{\sum_{k=0}^M b_kz^{-k}}{\sum_{k=0}^Na_kz^{-k}}$$

can be factored to the form:

$$\frac{b_0}{a_0} \frac{\prod_{k=1}^M (1-c_kz^{-1})}{\prod_{k=1}^N(1-d_kz^{-1})}$$

where each of the factors $(1-c_kz^{-1})$ contribute a zero at $z=c_k$ and pole at $z=0$. Similarly $(1-d_kz^{-1})$ contributes a zero at $z=0$ and a pole at $z=d_k$.

How is the above factorization done?
Is this the only representation?


texts: How to set threshold value for Hough Transform

texts:
I am new to Computer Vision, I am writing my own code for Hough Transform I have found the following algorithm:

Algorithm:

• The cell (i,j) corresponds to the square associated
with parameter values (θj, ρi).

• Initialize all cells with value 0.

• For each foreground point (xk,yk) in the thresholded
edge image
– Let θj equal all the possible θ-values

• Solve for ρ using ρ=x cos θj +ysin θj

• Round ρ to the closest cell value, ρq

• Increment A(i,q) if the θj results in ρq

• After this procedure, A(i,j)=P means that P points in
the xy-space lie on the line ρj=x cos θj +ysin θj

• Find line candiates where A(i,j) is above a suitable
threshold

I have written the code to increment the indices of hough accumulator array.

My Code:

%%
% Hough Transform

[rows,cols] = size(E); % size of edged image

rhomax = ceil (sqrt(rows^2 + cols^2)); % maximum value of rho

rhomax = rhomax * 2; % multuplied with 2 to cater both +rho and -rho values

accu = zeros(rhomax, 101); % initialization of accumulaor array with rhomax and 101 possible theta values

for i =1:1:rows % for all rows

    for j =1:1:cols % for all columns

        if(E(i,j) ~= 0) % if the pixel is ONE (part of the edge) then process it

            for theta = 85.0:0.1:95.0 % for the theta values 85.0 - 95.0

                % computing the theta index for hough accumulator
                theta_index = 0;
                for k =85.0:0.1:theta
                    theta_index = theta_index + 1;
                end
                %---

                rho = round (i*cos(theta) + j*sin(theta)); % computing the value of rho

                % computing the rho index for hough accumulator
                rho_index = 1;
                if(rho < 0)
                    rho_index = rho_index + abs(rho);
                elseif(rho > 0)
                    rho_index = rho_index + (rhomax/2) + rho;
                elseif(rho == 0)
                    rho_index = rho_index + (rhomax/2);
                end
                %---

                accu(rho_index,theta_index) = accu(rho_index,theta_index)+1; % incrementing the corresponding index in accumulator

            end

        end

    end

end

Now I want to proceed further with the threshold value. I will be very grateful if someone tell me how should I select my threshold value ?


texts: Does the Nyquist frequency of the Cochlear nerve impose the fundamental limit on human hearing?

texts:
The bandwidth of human hearing by empirical data is $20 \; Hz$ to $20 \; kHz$. A cochlear implant stimulates the auditory or acoustic or Cochlear  nerve directly so that the hearing can be improved in the case of stimulation mechanism upstream of the Cochlear nerve has degraded.
Let us assume that the ear mechanism has not degraded (such as in a young and healthy adult). The Cochlear implant can likely improve hearing, even in this case, by increasing the bandwidth by amplifying the effect of the ear drum vibration (sensor actuation). However the neurons connecting the Cochlear nerve to the hearing region of the brain have an upper limit on the sampling rate on the order of $1 \; kHz$.

Does the the Nyquist sampling theorem limit the superhuman hearing and sound localization capability made possible by a Cochlear implant?


texts: FIR design for arbitrary magnitude and phase: Why can I separate real and imaginary parts like this?

texts:
In this question about the design of a FIR filter with arbitrary magnitude and phase specifications user robert bristow-johnson suggested to split the desired complex frequency response $H$ into its real and imaginary parts $H_R$ and $H_I$, then use firpm or firls to design two linear phase filters for both, setting the ftypeoption to hilbert for the imaginary part. Apparently, the sum of the two resulting FIR filters would yield a filter that meets the specification.

I'm trying to understand why this works. Here's my attempt:
$$H_R = |H|\cos(\arg\{H\})$$
$$H_I = |H|\sin(\arg\{H\})$$
Now we design two FIR filters with firpm or firls, setting the ftypeoption to hilbert for the imaginary part. To my understanding, the frequency responses of the two resulting filters $\tilde h_R$ and $\tilde h_I$ will then look like this:
$$\tilde H_R = \Re\{H\}\ e^{jp(k)}$$
$$\tilde H_I = \Im\{H\}\ e^{jp(k)+j\frac\pi2}$$
Where $p(k)$ is the linear phase term and the second filter has the additional phase shift due to the hilbert option that makes it a type III or IV FIR filter. The sum of the coefficients $\tilde h[k] = \tilde h_R[k] + \tilde h_I[k]$ would then be a FIR with the frequency response
$$\tilde H = \tilde H_R + \tilde H_I$$
because adding the coefficients is basically the same as adding the filter output signals for the same filter input signal.

\begin{align}\tilde H =& \Re\{H\}\cos(p(\omega)) + \Im\{H\}\cos(p(\omega) + \frac\pi2) \\&+ j\left[\Re\{H\}\sin(p(\omega)) + \Im\{H\}\sin(p(\omega) + \frac\pi2)\right]\\
=& \Re\{H\}\cos(p(\omega)) - \Im\{H\}\sin(p(\omega)) \\&+ j\left[\Re\{H\}\sin(p(\omega)) - \Im\{H\}\cos(p(\omega))\right]
\end{align}

This looks nothing like the response $H$ that I actually want. And I don't see how this could work. Is the original approach wrong or am I?