texts: Matrix form of 2D-DFT for a vectorized image

texts: I want to apply 2d DFT to a N by N image. However, image is vectorized such that it is NxN by 1. How can i find the matrix form of 2d DFT such that resulting vector from multiplication is the vectorized NxN by 1 2d DFT of my N by N image?


texts: Other time-frequency-plane tiling than STFT, DWT, ConstantQ-Transform: multiresolution STFT?

texts:
It is known that

a) the STFT gives a rectangular tiling of the time-frequency plane

b) the Wavelet transform gives a non-linear tiling (better frequency resolution for low-frequencies, and better time-domain resolution for higher-frequencies)

c) Constant-Q transform (such as NonStationaryGaborTransform) have a logarithmic scale for frequency bins (instead of linear with STFT) and have a time-frequency tiling like this (y-axis logarithmic):

Question:
Is there an transform like this:

i.e. like a normal STFT, but for which the FFT-size would change for successive time frames.
Example: if there is a transient, the FFT-size is small (512) to keep good time-domain resolution, then a few time-frames later the signal is rather stationary, so the FFT-size is higher (8192) to have good frequency resolution.
It's like an adaptative / multiresolution STFT.
Last but not least: of course it's always possible to perform FFT with different window-length on successive time-frames, so the forward transform I'm describing is probably easy to do.
What I'm asking here is a transform which is invertible, i.e. in which we can compute an inverse transform to recover the initial signal.
Does this exist in Matlab or Python?


texts: Derive Maximum of zero phase rotated wavelet from rotated wavelet

texts:
I have one question concerning phase rotation of signals and I am not sure wether it works or not:

Let's assume we have a zero-phase wavelet x(t) i.e. Ricker wavelet.
Unfortunately, in our experiment, the wavelet y(t) we use is not zero phase anymore but rotated by an phase angle $\Psi$. Yet, we still can clearly identify the peak of the wavelet at x(T).

Question:
Is it possible to derive the position of the peak of the corresponding zero-phase wavelet $x(t)$ from $y(\tau)$?

My first steps are:
1) equation Ricker wavelet:
$$
x(t)=(1-2 \pi^2 f^2 t^2) e^{-\pi^2 f^2 t^2}
$$

2) equation rotated Ricker wavelet

$$
y(\tau)=(1-2 \pi^2 f^2 \tau^2) e^{(-\pi^2 f^2 \tau^2 - \Psi \tau)}
$$

My assumption was setting $y(t) = x(t)$ and solve for $t$ since $\tau,f,\Psi$ are known. Unfortunately the amplitude of the wavelet changes with the phase rotation.
How can I account for this?

please help!


texts: Analyzing an undersampled recording from an antenna with matlab?

texts:
I'm trying to analyze the power spectral density of a signal recorded from a Ham radio antenna. The signal is about 30 seconds long, recorded at a 48 kHz sample rate. The frequency is center frequency of 8 kHz recording, with the receiver centered on 5.275 MHz. The signal I'm currently looking at is supposed to be around 5.260 MHz, however I'm getting inconsistent results with the matlab script below:

[FileName,PathName] = uigetfile('*.wav','Select the wav file');
PathName = strcat(PathName, FileName);  %Adds the filename to the end of the pathnames

[x,Fs] = audioread(PathName);
[pxx,f] = pwelch(x,500,300,500,Fs);
plot(f,pxx)
xlabel('Hz'); ylabel('dB');

When I run this code with '5.260MHz.wav', it produces the following plot:

These frequency peaks all occur between 0-4 kHz, which is unexpected.
Since 5.260e6 mod 48e3 = 28e3, I would expect the signal to be at 48e3 - 28e3 = 20e3, or 20 kHz.

I created my own signal and then undersampled it as follows:

% signal properties
f0 = 5.260e6;           % frequency 5.260 MHz
sf = 8e3;               % bandwidth 8 kHz

% sample properties
Fs = 24e6;              % sampling frequency 24 MHz
N = 10e6;               % number of samples

% generate random frequency modulated sinusoidal signal
fi = smooth(randn(N, 1), 11);
fi = fi / std(fi) * sf + f0;
x = cos(2 * pi * cumsum(fi) / Fs);
Fs = Fs / 500;
x = x(1 :500: end);

filename = 'test_5.260MHz.wav';
audiowrite(filename,x,Fs);

When I plot this signal with the above script, it looks as expected, with a peak at 20 KHz:

Finally, I used the following script from the mathworks site to create another signal, which also graphed as expected:

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sample time
L = 1000;                     % Length of signal
t = (0:L-1)*T;                % Time vector
% Sum of a 50 Hz sinusoid and a 120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));     % Sinusoids plus noise

filename = 'test_50_120Hz.wav';
audiowrite(filename,y,Fs);

Because I'm getting consistent results with my generated files, I would like to believe that there is simply an error with the information concerning my antenna signal. However, I have several different antenna signals, and they are all graphing in the 0-4 kHz range, regardless of where the signal supposedly exists or the central frequency. This leads me to believe that I may be using the matlab code incorrectly, or that I may be misinterpreting how .wav files and radio signals work exactly.

Any input is appreciated. As I am still new to stackexchange, any suggestions on how to make the question clearer are also appreciated.


texts:
My master’s dissertation subject is to analyze and implement an algorithm that would sample (store) thermal images from camera with different frequency rates to optimize database size and processing power.
The thermal images are provided to computer with constant flow with packets structure consisting of header and matrix of 14 bits values per each ,,pixel''.
From publications that I found, most commonly compression techniques are used, but here I should develop dynamic sampling algorithm without althering data structrure.
Is anyone aware of any publication related to this subject? It seems that I cannot find any information on sampling patterns for raw 2d data streams, nor algorithms that would detect or calculate if the image differs enough from previous ones and is worth capturing.


texts: How to obtain impulse response from the differential equation of a system?

texts:
This year I'm having trouble with my Signals and Systems class. My major subject is Software Engineering and Electric and Electrical Engineering is my Minor. This question was my previous exam question but I coudn't write anything about it as a solution. If somebody could at least show a way to solve this problem it would be great.

Question:
$$\frac{d^2y(t)}{dt^2}+\frac{6dy(t)}{d(t)}+8y(t)=2x(t)$$
a) Find impulse response $h(t)$
b) Find output $y(t)$ if $x(t)=te^{-2t}u(t)$

I know that $h(t) = \frac {y(t)}{x(t)}$, and i m guessing that i can get that using z-transform but I don't know how to and if it is the correct way.

Thanks in advance...

EDIT:
Here is what I have done so far.

From here, to find the output I think I will use convolution.? y(t) = h(t)*x(t)

I used Laplace transform to find the inverse fourier transform of the function H(jw). What was I going to do if Laplace transform would not be suitable to situation?


texts: Rigorous proof of the convolution- Sum Representation of LTI Systems

texts:
We've started learning LTI systems and I have confusion regarding the convolution sum representation of a system.

My problem stems from insufficient mathematical formalism in signals and systems introductory textbooks.

We have the discrete-time unit impulse function (or kronecker delta) $\delta[n]$:

$$\delta[n] = \begin{cases}
 1 \qquad & \text{if } n = 0,\\
 0 \qquad & \text{if } n \neq 0 \\
\end{cases} $$

then $x[n]$ can be represented as:

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

I'm not sure how infinity is defined here (obviously as a limit, but I couldn't find a book where it's mentioned.). However, I understand the equation intuitively.

Then, as a 'direct' result from linearity and time invariance:

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

Where $h[n-k]$ is the impulse response.

Now, everything makes sense intuitively when I think of infinity as a "large enough index". But formally, infinity is a limit, and linearity is defined for a finite number of terms.

How can this be proven rigorously using a formal definition of infinite sum?


texts:
My efforts of solving this question are below.

I came to a conclusion that this system is causal, since:
$$
\begin{cases}
w[k]+5w[k-1]+6w[k-2]=x[k] \\
y[k]=w[k]+2w[k-1]+3w[k-2]+4w[k-3]
\end{cases}
$$
Where $w[n]$ is the output of the leftmost summation. After conducting Z transforms:
$$
\begin{cases}
\frac{W(z)}{X(z)} = \frac{1}{1+5z^{-1}+6z^{-2}} \\
\frac{Y(z)}{W(z)} = 1+2z^{-1}+3z^{-2}+4z^{-3}
\end{cases}
$$
$$\Rightarrow H(z) = \frac{Y(z)}{X(z)} =
\frac{Y(z)}{W(z)} \frac{W(z)}{X(z)} =
\frac{1+2z^{-1}+3z^{-2}+4z^{-3}}{1+5z^{-1}+6z^{-2}} =
\frac{z^{3}+2z^{2}+3z+4}{z^{3}+5z^{2}+6z} =
\frac{(z+1.651)(z+(0.175-j1.547))(z+(0.175+j1.547))}{z(z+2)(z+3)}
$$
According to that, I had concluded that this system has 3 zeroes and 3 poles, therefore the number of zeroes is less than or equal to the number of poles, and therefore the system is causal.
However, the answer key states the following:

Where have I made a mistake? Please, correct me if I'm wrong.


texts: Synthetic sounds to describe motion in 3d or 100d

texts:
Watch a butterfly flitting about, or an optimizer chugging along in 3d or 100d --
a sequence of points $X_1 \ X_2 \ X_3\  \ldots$

How could one generate synthetic sounds that convey

moving fast or slowly
moving steadily in one direction, or oscillating ?

The question has two parts (which should perhaps be asked separately):

how to filter a high-dimensional path to a 1d signal that conveys the kind of motion
how to convert that to pleasant sounds.

The goal would be to listen to sounds made by a running optimizer
as it moves along in 3d or 100d,
and glean some info on how it's moving.


texts: Library/strategy for rotate & scale image alignment/registration

texts:
I have two raster images of maps, A and B.  I know that B contains A but I do not know the scale or rotation.  There should not be any other distortions between the image but there may be noise.  Please could anyone recommend an approach or software library which would provide me with the scale and rotation transformations to align A onto B, thus allowing me to transform arbitrary points between the coordinate spaces?

Left: A: source, B: destination.

Right: Location of A on B

Given the high contrast, plane sharing and 'edgyness' of the images I would hope it would be fairly trivial compared to most image analysis tasks.  There are similar questions on here but they seem to be doing something more complex.


texts:
I'm trying to decimate given picture

Vertically.
Horizontally.
Horizontally + Vertically.

I know that I should use low-pass filter* and then downsample** it not to see any aliasing in image.
To this my code is below:

Horizontally filter:

filter_hor    = [0.25,0.5,0.25];
filtered_hor  = imfilter(big_shirt,filter_hor);%% big_shirt is shirt image given above
decimated_hor = filtered_hor(1:2:end,:);%%decimated_hor = downsample(filtered_hor,2);
imshow(decimated_hor)
title('Horizontally filtered')
axis on

Vertically filter:

filter_ver    = [0.25;0.5;0.25];
filtered_ver  = imfilter(big_shirt,filter_ver);
decimated_ver = filtered_ver(:,1:2:end);
imshow(decimated_ver)
title('Vertically filtered')
axis on

Horizontally + Vertically filter:

filtered_2D  = imfilter(decimated_hor,filter_ver);
decimated_2D = filtered_2D(:,1:2:end);
imshow(decimated_2D)
title('2D filtered')

What I get is :

*My low pass filter is: [0.25,0.5,0.25]

**By a factor of 2

My expectation is:

When I use horizontal filter, I just get rid of aliasing on lines of shirt collar.

When I use vertical filter, I just get rid of aliasing in vertically lines.

Vertical+horizontal will be nice picture like what I get.

However in first two cases are not satisfied. What can be the problem. My resolution on computer is 1366X768 and image size is 800X600. When I use imshow, I get:

Warning: Image is too big to fit on screen; displaying at 67%

Can be problem related to this or am I doing something wrong?


texts: how to robustly estimate low and up envelope of signal with trend, few level constant steps and noise

texts:
I am looking for robust estimation method of low and up envelope of the signal consisting from smooth trend component, constant steps between few fixed levels and additive noise (+ outliers of course). This question raised during my current research work on real life signal processing.

Typical example of signal is produced by following MATLAB script:

  %% signal definition
  % number of samples
  Ns = 10000;
  % sampling period [secs]
  Ts = 1;
  time = (1:Ns)*Ts;
  % trend component
  a = 2;
  b = 0;
  T = 1e4;
  slope = time/Ts * 0.0005;
  trend = a * sin(2*pi*time/T) + b + slope;
  % steps component (4 constant levels)
  step = [zeros(1,1000),linspace(0,1,6),1*ones(1,500),linspace(1,0,6),zeros(1,500),linspace(0,-2,10), -2*ones(1,1200),linspace(-2,1,15),1*ones(1,3000),linspace(1,0,6), zeros(1,3000),linspace(0,-1,6), -1*ones(1,751)];
  % noise component (normal noise)
  noise = 0.2*randn(1,Ns);
  % noise = 0.5*(rand(1,Ns)-0.5);
  %
  %% show signals component
  close all
  figure
  plot(time,trend,'r-')
  hold on
  plot(time,trend+step,'g-')
  plot(time,trend - 2, 'b--')
  plot(time,trend + 1, 'k--')
  plot(time,trend+step+noise,'mo')
  legend('trend','trend+steps','lowenvelope', 'upenvelope','trend+steps+noise')
  title('smooth trend signal with constant steps between 4 levels and noise')
  xlabel('time [sec]')
  ylabel('value [-]')
  hold off

See the following image

The separate signal components are unknown! Steps are always constant and between  small number of fixed levels (typically < 4 or 5), so the estimated envelopes should be parallel to the trend signal. Noise is approximated by normal distribution with sigma ~0.1

Any idea how to solve this surprisingly difficult problem? Any relevant references or matlab code links?

Add notes:

In principle I am looking for any robust estimation of trend signal component which is nearly parallel to real trend signal.


texts:
To estimate the doppler shift in a multiuser communication, I use preambles (PN codes). When I send the preambles with the datas (BPSK symbol), the estimation for the doppler shift are not correct. But with just the preambles and no datas the estimations are correct.

Where is the problem?


texts:
I am currently trying to understand the Chirp Z Transform, and compute an FFT for a vector that has size which is not necessarily a power of 2. The key step in the Chirp Z transform to recognize that the boils down to computing a circular convolution of two vectors of length $n$.
My question is: How does one compute the circular convolution of two vectors say ($x[1],x[2],...,x[n]$) and ($y[1],y[2],...,y[n]$) using DFTs only sizes that are a power of 2? I know that by the convolution theorem, I can find the circular convolutions of two vectors applying the DFT (and it's inverse) a few times, but I only have power of 2 sized DFTs at my disposal to compute this for arbitrary $n$. I considered padding the vectors to increase their size to the nearest power of 2, but it seems as if one must carefully construct such a padding to get the right result (since we need to compute modulo $n$ and not $2^k$).
For more context: I was going through Dilip Sarwate's answer on this topic at this link: https://math.stackexchange.com/questions/104148/chirp-transform-and-convolution[on math stackexchange]1 and was unable to figure out how one can compute a circular convolution efficiently.


texts: Overlapping in real time fourier transform?

texts:
I have an algorithm and I need to record audio and perform short time Fourier transform to obtain which frequency is the most common. I am using a Hanning window to try and reduce spectral leakage as much as possible.  My problem is how should I overlap the data because since it is in real time I will no have the next part of the data.  Should I perform the Fourier transform in delay so that I can overlap the next segment.
So basically I am computing the the Fourier Transform of the red part, I do not have the blue part yet because it is still in the future.  How can this problem be solved?  Should I wait for the next segment to be recorded and perform the Fourier Transform of the past data so effectively having a delay?


texts: Save computational burden by calculating the output of decimate filter at output rate?

texts:
I read in a book that the output of a decimate filter can be calculated at output sampling rate. In this way the number of filtering operations can be greatly reduced. But I am wondering how to implement that in practice.
The built-in filter functions provided by Matlab and a lot of realtime DSP libraries are optimized to have best computation and storage performance, but they don't provide the capability of calculating the filter output only at output sampling rate. Is it possible to write a code, in a for-loop style, that only calculates the decimated data, and outperforms the performance of built-in functions? If not, what modifications can be done to the built-in functions to reduce the computation of decimation filter, especially in realtime implementations?
Thank you


texts: How do I add more resolution / taps to an analog -> digital filter?

texts:
(I am trying to create an IIR audio filter that adds reverb to an initial sample)

Say I designed an analog filter to model acoustic attenuation based on the following mathematical model:

$$
I = I_0 e^{pt},
$$

Where $p$ is some constant such that $-1 \leq p < 0$.

In the laplace domain, this is simply

$$
\mathfrak{L}\{I\} = \dfrac{1}{s-p}
$$

Using the impulse invariant transform, I get the digital equivalent

$$
Z\{I\} = \dfrac{z}{z-e^{pT_s}},
$$

Which inherently has one pole.

How do I add more resolution to this filter? Implementing this as a digital filter makes my audio sample sound like garbage. Would adding arbitrary taps keeps it from modeling the original analog decay?

Using experimental attenuation data and frequency sampling for an FIR equivalent, I can easily obtain 50000 poles and create a very clear reverb-adding filter.

Both IIR and FIR methods were tested in MATLAB, with the FIR using MATLAB's built-in church impulse response.

(Sorry if the question is unclear, I'm a bit new at this.)


texts:
I gotta wrap my head around to design CIC compensation filter

I'm studying by referring these materials:

Altera, "Understanding CIC compensation filters
Hardware Efficient FIR Compensation Filter
For Delta Sigma Modulator Analog to Digital Converters.
Circuits and Systems, 2005. 48th Midwest Symposium on
Saiyu Ren ; Dept. of Electr. Eng., Wright State Univ., Dayton, OH ;    Siferd, R. ;    Blumgold, R. ;    Ewing, R.

first Altera just say "inverse sinc function"
they wouldn't let me know how to realize the function except for using matlab function fir2 and how to get the coefficients of taps

so I found another reference (2) which I mentioned above.

they explain how to decide the coefficients of taps. but I couldn't figure out some equations.

they say

H(z)=a0*z^0+a1*z^-1 + ... + an*z^-n  => H(f)=h(0)+2*sigma(from 1 to (n-1)/2) h(k)*cos(2k*pi*F) where F=f/fs

Here is my first question. How is H(z) converted to H(f). If I substitute z to exp(iwT) it doesn't make sense.

second question I though a0=an=h(1). Is it correct?

finally i used matlab to calculate the coefficients of taps and i got the same coefficients in comparision with a result of the paper, but a simulation gave me a totally different and wrong data

h= [-0.70  2.09 -1.76  0.74 -1.26 0.40  0.34  0.26 0.74
    -0.53  0.24 -0.74  0.04 -0.39 0.05 -0.02  0.28 0.14
     0.41 -0.02  0.36 -0.02  0.41 0.14  0.28 -0.02 0.05 -
     0.38  0.04 -0.74  0.24 -0.53 0.74  0.26  0.34 0.40 -
     1.26  0.74 -1.76  2.09 -0.70] 41 taps

please let me know where I have to modify...

If I violate The copyright issue (because of the content of paper)
please tell me and I will edit or delete this.

Here's my matlab code and plot

clc; clear;
OSR=8
fb=125e6/2

fs=2*fb*OSR
f=[0 14 28 38 45.5 49.5 100 140 170 200 225 250 275  300 325 350 375 400 430 460 490]*10^6;

F=f/fs;
N=41;
L=[1: 1 : (N-1)/2];
S=max(L)

for i=1:S+1
   A(i,1)=1;
end

for i=1:S+1
    for j=2:S+1
    A(i,j)=2*cos(2*pi*(j-1)*F(i));
    end
end

H_sinc=((sin(OSR*pi.*F)./(OSR*sin(pi.*F)))).^2;

H_FIR=1./H_sinc;

H_IFIR(1:6)=H_FIR(1:6);

H_IFIR(1)=1;

H_IFIR(7)=0.1;

H_IFIR(8:21)=0;

H_IFIR=H_IFIR';

A_inv=inv(A);

H=inv(A)*H_IFIR;

H=H';

for i=1:S+1
  H_sol(i)=H(S+2-i);
end

for i=1:20
    H_sol(N-i+1)=H_sol(i);
end
f=[0:100:fs];
F=f/fs;

H_fir1=H_sol(21)+2*H_sol(1)*cos(2*pi.*F) ...
           +2*H_sol(2)*cos(2*pi*2.*F) ...
           +2*H_sol(3)*cos(2*pi*3.*F) ...
           +2*H_sol(4)*cos(2*pi*4.*F) ...
           +2*H_sol(5)*cos(2*pi*5.*F) ...
           +2*H_sol(6)*cos(2*pi*6.*F) ...
           +2*H_sol(7)*cos(2*pi*7.*F) ...
           +2*H_sol(8)*cos(2*pi*8.*F) ...
           +2*H_sol(9)*cos(2*pi*9.*F) ...
           +2*H_sol(10)*cos(2*pi*10.*F) ...
           +2*H_sol(11)*cos(2*pi*11.*F) ...
           +2*H_sol(12)*cos(2*pi*12.*F) ...
           +2*H_sol(13)*cos(2*pi*13.*F) ...
           +2*H_sol(14)*cos(2*pi*14.*F) ...
           +2*H_sol(15)*cos(2*pi*15.*F) ...
           +2*H_sol(16)*cos(2*pi*16.*F) ...
           +2*H_sol(17)*cos(2*pi*17.*F) ...
           +2*H_sol(18)*cos(2*pi*18.*F) ...
           +2*H_sol(19)*cos(2*pi*19.*F) ...
           +2*H_sol(20)*cos(2*pi*20.*F);

figure(1), semilogx(f,db(H_fir1))
% hold on
grid on

please help me to solve this problem...


texts: how important is the pre-emphasis stage in MFCC?

texts:
I wrote a simple MFCC routine, using all the steps usually suggested in the literature with the exception of the pre-emphasis stage. How important is this stage?
 Also, some authors suggest that using rectangular windows instead of triangular doesn't degrade MFCC performance (as long as the windows are overlapping) so any ideas regarding this particular subject are welcome.


texts: Power allocation for frequency selective channel [waterfilling ]

texts:
I am studying channel capacity with channel state information at a transmitter, where the water-filling algorithm is discussed. The algorithm should allocate more power to some subcarriers with larger SNR for maximizing the channel capacity in systems, as I understand. This algorithm discusses in case flat MIMO channel. For me interesting how the water-filling algorithm works in case frequency selective channel.

Is there water-filling algorithm a difference when it is computed for flat or frequency selective channel? if yes, what has to be changed in the simulation? Water_level (u)?

Matlab implementation:

    function s = WaterFilling(lemda, SNR)

        lemda = sort(lemda, 'descend');
        n = length(lemda);
        s = zeros(1, n);
        loop = 1;

        i = 1;
        while(loop == 1)
            loop = 0;
            sum = 0;
            for k = 1:n-i+1
                sum = sum +1/(SNR*lemda(k));
            end
            u = 1/(n-i+1)*(1+sum);  % water-level
            for k = 1:n-i+1
                s(k) = u-1/(SNR*lemda(k)); % power allocation, gamma
                if s(k)<0
                    loop = 1;
                    s(k) = 0;
                end
            end

            if loop == 1
                i = i+1;
            end
    end


texts: I have seen this be mentioned a couple of times in some books I read, so I want to make sure. Is the complex envelope simply the summation of the real and quadrature components of a signal, whereby the absolute value is the (real) envelope? I have read this wiki page but I am not sure I completely understand. Is the complex envelope simply the real and imaginary parts of a passband signal combined? Thank you.


texts: What is the difference between ISI and MAI in MIMO system?

texts:
ISI ： inter-symbol interference

MAI ： Multiple Access Interference

I know that multi-path would caused delay spread, and when signal symbol period is less than delay spreaad, ISI would happened.
In MIMO system, can I say that owning to the differenct location of antenna element that caused different path length, which make different delay and caused MAI ?
I also want to ask :
When would MAI happend?


texts: Coupling two resonant bandpass filters? (To simulate guitar/piano string effects)

texts:
A major part of the sound of a guitar or piano string comes from the coupling of the horizontal and vertical string vibrations. This concept is described in terms of transfer functions here.

This article describes how the interactions between the horizontal and vertical vibrations can be modeled as two strings coupled by an ideal transformer:

The transformation ratio is estimated for vertical > horizontal transfer to be 5:1:

  The transformation ratio n:1 of the vertical  to horizontal string velocities at the bridge, which determines the coupling degree, is assumed to be n=5, the value for a grand piano measured by Mori et al.[9].

Mori et al. states here:

  The  amplitude  of  the  horizontal  reflected  wave  is  approximately  1/10  of  the  vertical  input  amplitude.  From  this  result,  it  can  be  considered  that  the  transform  ratio  of  reflection  is  5:1,  because   the   observed   wave   consists   of   the   direct   component  and  the  component  reflected  by  the  bridge;  The ratio of input impulse, the vertical reflection and the horizontal  reflection  is  5:4:1.

I am modeling strings using arrays of resonant bandpass filters. I intend to have one array for the horizontal and one array for the vertical vibrations in this manner and couple them as described.

I am wondering in a very basic sense - how do you couple bandpasses like this? In the simplest case, let's just imagine one bandpass for the horizontal and one bandpass for the vertical. I excite them both. Would I then just feed a small percent of the output of one (delayed by one sample) into the input of the other and vice versa? Would I attenuate each bandpass a bit to compensate for the loss of energy from coupling so they don't become infinitely sustaining?

Eg. Here is what I am thinking in C++:

bandpass1.setQ(1/((1/bandpass1Q) + (1/couplingQ)));
bandpass2.setQ(1/((1/bandpass2Q) + (1/couplingQ)));

double bandpass1Output = bandpass1.nextSample(exciter + bandpass2Output_z1);
double bandpass2Output = bandpass2.nextSample(exciter + bandpass1Output_z1);

double bandpass1Output_z1 = bandpass1Output * 0.1;
double bandpass2Output_z1 = bandpass2Output * 0.1;

Does that make sense? In principle, I am adding a bit of damping (couplingQ) to compensate for the energy lost from coupling. Then I am taking some fraction of the prior sample's output from each bandpass and feeding it into the other bandpass.

Is this the correct principle? If so, what percent of each would I be feeding from one to the other at each sample to simulate a 5:1 transformation ratio as described? How would I calculate the couplingQ from this as well?

Any guidance would be appreciated. Thanks.


texts:
I'm coding a synthesis external for MaxMSP in C.
It's based on Peter Blasser's bounds and bounce concept ( http://petermopar.blogspot.co.uk/2014_04_01_archive.html )- basically triangle waves with constant cursor rate (rather than frequency) and modulated bounds - eg instead of moving between -1 and 1, the wave moves between lo and hi, with lo and hi modulated.

This morning I've been thinking of adding waveshape morphing between triangle and sine (or at least sinusoidal since the moving bounds will deform the sine).

I can think of a fairly simple way to do this:

// Assuming classic triangle moving between bounds -1 to 1
// and defining position of triangle cursor @ n = pos
// shape_amt between circa 0.1 (Tri) and 1.0 (sine)
out = sin(pos * shape_amt * 0.5* pi ) / sin(shape_amt *0.5 * pi())

ie - basic waveshaping technique - scale the portion of the lookup table that's being scanned (in this case some segment of half a sine wave) to vary the shape.

For my specific case all I'll need to do from there is build in scaling from relative position between to bounds to relative position between -1 and 1.

My question really is- is this a horribly inefficient way of achieving what I'm after? I'm new to DSP coding and to coding in C so I'm sure the way that's most obvious to me isn't the best.

The obvious performance improvement that occurs to me is to create a wavetable of the sin function and look up rather than calculate those values. Obviously that way I lose resolution on the "straighter" portions of the wave, but either a high res lookup or decent interpolation feels like it ought to deal with that fairly well. Mainly I'm considering this for LFO applications so quality needn't be utterly stellar (though it'll be nice nice if it's useable into the audio range without too much fuzzy aliasing).

But beyond that, is there just a completely different approach that would be more efficient?
And if this is the best approach, are there any things I should think of when pursuing it?

Edit - Thanks for great answers below. I've ended up doing something fairly simple (and possibly crude) for the moment, which suits my purposes pretty well. I'll definitely be investigating both suggestions as improvements to this.

#define sign(a) ( ( (a) < 0 )  ?  -1   : ( (a) > 0 ) )
#define LKTBL_LNGTH 2048

void setup_lktables (t_object* x)
{   // create lookup for 1/4 sine
    int i;
    for(i=0, i< LKTBL_LNGTH, i++){
        x->sin[i] = sin(pi * i * 0.5 / LKTBL_LNGTH) ;
    }
}

double lookup_shape (t_double pos, t_double shape, t_double lo, t_double hi)
// acting on triangle wave (phase + (2*gradient*f/sr) per sample, wraps & changes dir  at <lo and >hi)
// shape comes in as restricted to (-1...-0.05, 0.05 ...1), defines the portion of lookup to use

{
    t_double midpoint, halfwidth, ph, fracph, shaped;
    t_int maxph, intph, sign;
    // get relative position between bounds for waveshaping lookup
    midpoint = lo + 0.5f * (hi - lo);
    halfwidth = midpoint - lo;
    // prepare phase values for lookups
    maxph = shape * LKTBL_LNGTH - 1;
    ph = (pos - midpoint) * maxph /  halfwidth;
    sign = sign(ph);
    ph = ph * sign;
    intph = (int)ph;
    fracph = ph - intph;
    // lookup, scale & lerp
    shaped = sign * (x->sin[intph] * (1.f - fracph) + x->sin[intph+1] * fracph) / x->sin[maxph];

    //now return  waveshaping output scaled to actual bounds
    return midpoint + shaped * halfwidth
}


texts: Signal response in Frequency Hz from angular frequency

texts:
I have a filter whose respose I have in angular frequency.

A0 = [1.0000000,0.3281250,0.0000000];
B0 = [0.5859375,-0.0625000,0.0000000];
[z,p,k]=tf2zp(B0,A0);
figure;
zplane(z,p);
w = 0:(pi/500):pi;
H=freqz(B0,A0,w)
figure;
plot(w,20*log10(abs(H)/max(abs(H))))

1) How can I get the corresponding response in Hertz? Is it:

plot(w/(2*pi),20*log10(abs(H)/max(abs(H))))

2) Also, how do I get the response of the above filter, if I apply it to a sinusoidal signal at 20MHz, with a sampling freq of Fs=300MHz


texts: How does google map put together satellite pictures?

texts: Just out of curiousity, does anyone know how does google map put together all those massive quantity of satellite pictures with with high accuracy? What advanced computer science technology has it employed?


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

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

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

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


texts:
I'm currently reading articles about the different type of beamformers. This spatial filtering is gonna be used for acoustic purposes, to focus the beam on the person we need to hear talking in a crowded environment.

I think that there is a misconception on some of the things I taught I understood.

So, the way I saw it, I should decide on a criterion to optimize (MVDR or MSNR for example), and once I had determined the optimal weight, I should use an estimator algorithm (like LMS/RLS) to converge to this solution.
More explicitely, the steps I had in mine were:

$\mathbf{w}_{\text{opt}} = \frac{\mathbf{R}^{-1}_x \mathbf{d}}{\mathbf{d}^H\mathbf{R}_x\mathbf{d}}$, $\mathbf{R}_x$ being the covariance matrix, and $\mathbf{d}$ the steering vector
$\mathbf{y}_{\text{desired}} = \mathbf{w}_{\text{opt}}\mathbf{x}$, $\mathbf{x}$ being the vector of our different entry
Use LMS with $\mathbf{y}_{\text{desired}}$ as our reference to estimate the weights of our beamformer

However, this article from ResearchGate (https://tinyurl.com/uoexq4c) talks about MVDR beamformer and LMS beamformer as two different things.

From this point, I don't understand how neither of them works. If I don't use the MVDR optimum solution, I can't have a reference signal for my LMS, and I can't find any way to minimize the MVDR criterion without using the optimal solution(which might be complex to use when i'll implement my algorithm on an STM32).

Therefore, I don't really understand how to use neither of these two methods separately.

How should I approach the optimal weights of my MVDR?
What is my reference signal for my LMS if i don't use MVDR optimal weights?

Thank you for your responses.


texts: What does the word Entropy refer to in Entropy encoding

texts: I have come know of Huffman and Arithmatic encoding that try to curtail redundant data and thus compress data. Why are they called Entropy encoding? What does that mean? Why not just call them compression encoding?


texts: Absolute and Square Summability & Integration of sinc function

texts:
Is the sinc function both absolutely summable (L1 norm for Continuous time signals and l1 norm for Discrete time signals) and square summable (L2 norm for Continuous time signals and l2 norm for Discrete time signals) ?

Can anyone show the integration and summation calculations that calculates above norms ?


texts: Carrier Frequency Offset Correction Limitations

texts:
I know that carrier frequency offset is typically detected and corrected in OFDM systems. I was wondering if there is a way to determine the theoretical limit to how much offset can be successfully "compensated for" and still utilize the data successfully (specific to OFDM, if that matters). I'm thinking in the simplest terms where the transmitter and receiver are trying to transmit and receive at the exact same rate (for example, the receiver isn't simply oversampling). I know it has to do with sampling rate and the overall bandwidth of the OFDM packet...

edit: maybe I should have specified, but I'm referring to CFO correction done in the digital domain, post ADC (think DSP) instead of a voltage-controller oscillator type of control loop.


texts:
I was working on wavelet signal decomposition and got confused with the downsampling part of Discrete Wavelet Transform algorithm.

If we consider a signal $\{a_0[0], ... ,a_0[N-1]\}$ of length N and use formula
$$
a_1[k] = \sum_m h[m-2k]a_0[m]
$$
then for two low pass filters
$$\{h[0]=h_0,~h[1]=h_1, ~h[2]=h_2, ~h[3]=h_3\}$$
and
$$\{g[-2]=h_0,~g[-1]=h_1, ~g[0]=h_2, ~g[1]=h_3\}$$
(which are the same thing except the second filter is time-shifted version of the first one) we will get different results.

Let me clarify. After convolution of $a_0$ and $h$ we will get $N+3$ samples $\tilde a_1[k]$, $k=-1,...,N+1$ with the first non zero sample being $\tilde a_1[-1]$. So, after downsampling we will result in the following sequence
$$
\tilde a_1[0], \tilde a_1[2], \tilde a_1[4], ..., \tilde a_1[N_1]
$$
where $N_1$ is the integer part of $(N+1)/2$, with total number of samples equal to $N_1 + 1$.

Now if we convolve $g$ and $a_0$ again we will get $N+3$ samples, but this time $\hat a_1[k]$, $k=0,...,N+2$. And after downsampling we will get
$$
\hat a_1[0], \hat a_1[2], \hat a_1[4], ..., \hat a_1[N_2]
$$
where $N_2$ is the integer part of $(N+2)/2$, with total number of samples equal to $N_2 + 1$.

So as a result, if number of samples $N$ of input signal $a_0$ is even, then number of approximation coefficients for $h$ and $g$ is different (integer part of $(N+1)/2$ is less than integer part of $(N+2)/2$). More than that, second sequence consists of exactly those samples which are excluded from the first one and vice versa.

So my questions are:

1) Is this correct or do I have a mistake somewhere?

2) If this is correct, then why implementations of dwt algorithm ignore this fact? (I've tested dwt from PyWavelets package)


texts: Down sample incoming stream of data to plot a dynamic graph

texts:
I am working on a module where I need to plot ECG streaming data on a graph. This graph has configurable speed i.e. plot 1 heartbeat data sample within say 50mm,25mm, etc. For this, I am required to down sample based on the screen resolution and plot the same.

I am using Qt QML for plotting. What is the exact algorithm to achieve this? I don't have a clue on what exactly I have to google


texts:
I am trying to understand how to interpret a spectrogram image plot: (1) what does the yellow band and its corresponding color code with a high positive value indicate? There are so many yellow bands in comparison to blue. The legend shows that the yellow bands have higher positive values, so is that higher frequency?
(2) I used the following code to generate the spectrogram. The X axis contains 50 data points because the method splits the signal into segments of length 50. I arbitrary used 10 into the spectrogram function. Is there a specific rule as to how to determine the number of points?
x(1) = rand;
for n = 1:500
    x(n+1) = 4*x(n)*(1-x(n));
end

specto=spectrogram(x,10,0,256); %using 256 fft points with no overlap

logspecto=(20*log10(abs(specto)));


texts:
I am trying to characterize a power amplifier fitted to transmitter board.
I applied a sinewave at the required carrier frequency and I collected the input samples and the output samples by using a digital scope. I have noticed a strange behaviour which I can not explain.

when the input sample is zero the output is not zero but a significant value.

Should I expect when the input sample is zero, the output sample from the power amplifier should be zero ?

May the input and the output samples are not in sync.
I tried to sync them, but still the zero input samples do not correspond to very small signal or zero sample.

When I plotted the output/input  I received this plot.

Is this an expected plot ?

How I can get a smooth plot for the ratio of the output over the input ?

I appreciate any comments.


texts: Downsampling fft results to a graphics display

texts:
I'm new to dsp.  I've got a complex signal that was sampled at 192kHz then an fft of 16384 bins created.  This now needs to be plotted on a graph for a real time spectrum display.  The problem is if you want to plot the whole spectrum consisting of a 192kHz span, you will be plotting 16384 points on a display that many times will be 1600 pixels wide for example.  The result is you get an ugly mess on the screen.  So far I have been downsampling the signal by just plotting every 8th point. The problem is though, when I look at other spectrum displays that also perform 16384 fft's, they show much more data.  I have effectively dropped some on the peaks that the other software doesn't.
So the question is, how does one do the downsampling properly?
Thanks


texts:
Introductory resources to mathematical properties of "dynamically recalculable" filters in audio/musical equalizers?

Particularly, I want to understand

What mathematical features (e.g. monotonicity of coefficients) of filters are suited for real-time modulated filters?
How to analyze filters for the required features?

But I want to also understand things such as "filter instability" and the ways they are implemented in real-time (does one apply them in buffers or what?)?

Dynamic recalculation: That the filter's cutoff and resonance can be recalculated on the fly. Or that the parameters can be adjusted while the filter is running.


texts: Is this the correct solution for finding the samples of phase delay at a given frequency of an IIR BiQuad filter?

texts:
I want to calculate the number of samples of phase delay at a given frequency of an IIR BiQuad Filter. I believe I have solved this, but I would like to confirm I have done it correctly. My work is as follows.

Background:

The transfer function of a BiQuad filter is:

$$
H(z) = g \frac{b_0 + b_1z^{-1} + b_2z^{-2}}{a_0 + a_1z^{-1} + a_2z^{-2}}
$$
ω is normalized by the sample rate (samples/sec) with units of radians/sample. The normalized Nyquist frequency is π radians/sample and the normalized sample rate is 2π radians/sample.

$$
\omega = \frac{2 \pi f}{f_s}
$$
The frequency response of a system can be calculated by letting $z = e^{j\omega}$, where ω is a normalized frequency in the range (−π,π):

$$
H(e^{j\omega}) = g\frac{b_0 + b_1e^{-j\omega} + b_2e^{-j2\omega}}{a_0 + a_1e^{-j\omega} + a_2e^{-j2\omega}}
$$

The value of $H(e^{j\omega})$ will be complex. The phase response at the frequency ω is just the phase angle of the resulting complex number. The magnitude response at the same frequency is likewise equal to the magnitude of the number.

Euler's Formula for solving complex exponents is:
$$e^{ix}=cos(x)+isin(x)$$

Solution:
$$
H(e^{j\omega}) = g \frac{b_0 + b_1e^{-j\omega} + b_2e^{-j2\omega}}{a_0 + a_1e^{-j\omega} + a_2e^{-j2\omega}}
$$
$$
H(e^{j\omega}) = g \frac{b_0 + b_1(cos(-ω) + jsin(-ω)) + b_2(cos(-2ω) + jsin(-2ω))}{a_0 + a_1(cos(-ω) + jsin(-ω)) + a_2(cos(-2ω) + jsin(-2ω))}
$$
$$
H(e^{j\omega}) = g \frac{b_0 + b_1cos(ω) + b_1jsin(-ω) + b_2cos(2ω) + b_2jsin(-2ω)}{a_0 + a_1cos(ω) + a_1jsin(-ω)) + a_2cos(2ω) + a_2jsin(-2ω)}
$$
$$
H(e^{j\omega}) = g \frac{[b_0 + b_1cos(ω) +b_2cos(2ω)] + j[b_1sin(-ω) +  b_2sin(-2ω)]}{[a_0 + a_1cos(ω) + a_2cos(2ω)] + j[a_1sin(-ω))  + a_2sin(-2ω)]}
$$
$$
H(e^{j\omega}) = \frac{[g(b_0 + b_1cos(ω) +b_2cos(2ω))] + j[g(b_1sin(-ω) +  b_2sin(-2ω))]}{[a_0 + a_1cos(ω) + a_2cos(2ω)] + j[a_1sin(-ω))  + a_2sin(-2ω)]}
$$

Simplify:

$$m = g(b_0 + b_1cos(ω) +b_2cos(2ω))$$

$$n = g(b_1sin(-ω) +  b_2sin(-2ω))$$

$$p = a_0 + a_1cos(ω) + a_2cos(2ω)$$

$$q = a_1sin(-ω))  + a_2sin(-2ω)$$

$$H(e^{j\omega}) =\frac{m+in}{p+iq}$$

Find the Angle (delay in samples?):

For a given complex equation of:

$z= x+iy$

Where φ is the angle between the x axis and the vector z measured counterclockwise in radians:

$φ = atan2(y, x)$

Additionally, the phase of the ratio of two complex numbers is the difference of the phases of the numerator and denominator (numerator - denominator), as noted by Matt L in his answer below:

$$φ(ω) = atan2(n,m) - atan2(q,p)$$

To get the final phase delay in samples, we must take the negative of the phase we have solved and then divide that by the normalized frequency:

$$\tau_p(\omega)=-\frac{φ(\omega)}{\omega}$$

Therefore the final equation for phase delay for any normalized angular frequency ω would be:

$$\tau_p(\omega)= \frac{-atan2(n,m) + atan2(q,p)}{\omega}$$

(Updated from Matt L's feedback.)


texts: how to find the inverse response of a system

texts:
For the system described by the differential equation below find its inverse zero-state unit step response $$\dfrac{d^2y(t)}{dt^2}-2\dfrac{dy(t)}{dt}-8y(t)=\dfrac{d^2x(t)}{dt^2}-2\dfrac{dx(t)}{dt}-3x(t)$$
so that the cascade system of its inverse system and this system is the identity system.
I am not sure how to tackle this type of problem. My first thought is to switch $x$ and $y$ and find the response of that system using laplace transform, which I presume would be the inverse system. Is this reasoning correct? if not, how does one tackle this problem? I'm not asking for a solution (although one is welcome), but for a general approach to solve problems like this. Thanks!


texts:
Across the channel bandwidth there shoudd be 100 RBs (central 1200 sub carriers in a 2048 bin). Max power normalized to 0dB, the first OOB lobe should be -34dB. This 1st lobe on either side is at 5 MHz across the x axis from the symbol.

I am trying to get hold of filtering but I still have rookie doubts.

I need to apply a filter to 6 symbols of a time domain OFDM signal.
Each symbol is a 2048 FFT/IFFT bin having 1200 central subcarriers and has CP of 512 samples added (therefore total 2560 samples).

I am representing the PSD w.r.t FFT points - so 2048 samples?

While designing filters in matlab using fir1:

Do I select the cut off freq as 1200/2048 or is that wrong?
Do I filter all the symbols one by one in a loop without zero padding after every symbol gets filtered (for flushing the tail)?
Do I clone the same coeff's and filter in the receiver side?


texts:
I'm trying to solve a problem on convolution from Alan V.Oppenheim:
Find the convolution output $y[n]$ for the following signals:
$$x[n]= u[n]\quad\text{and}\quad h[n]=a^{n}u[-n-1], \ a>1 $$
I started the evaluation:
$$y[n]=\sum_{k=-\infty}^{+\infty} u[k]a^{n-k}u[-n+k-1]$$
considering that $u[k]=1$ for $k\ge0$ and $u[-n+k-1]=1$ for $k\ge n+1$
which I evaluated to
$$y[n]=a^{n}\sum_{k=m}^{+\infty} a^k$$ where $m=n+1$ could be $<0$ or $>0$
and I tried to evaluate for $m>0$ which is the same as $n>-1$, which evaluated as:
\begin{align}
y[n]&=a^{n}\sum_{k=m}^{+\infty} a^m\\
&=a^{n}\left[\left(\sum_{k=0}^{+\infty} a^m\right)-\left(\sum_{k=0}^{m-1} a^m\right)\right]\\
&=a^{n}\left[\left(\frac{1}{1-a^{-1}}\right)-\left(\frac{1-a^{-m}}{1-a^{-1}}\right)\right]\\
&=\frac{a^{n-m}}{1-a^{-1}}\\
&=\frac{a^{-1}}{1-a^{-1}}\tag{since $n-m=-1$}
\end{align}
but when I evaluated for $m<0$  which is $n\le -1$ I am facing a problem:
$$y[n]= a^{n}\left[\left(\sum_{k=m}^{-1} a^m\right)+\left(\sum_{k=0}^{+\infty} a^m\right)\right]$$
How do I evaluate the first summation? I mean am I to consider $k=-m$ since $m<0$?


texts: Why do linear systems show sinusoidal fidelity?

texts: I am looking for a proof for sinusoidal fidelity. In DSP we study a lot about linear systems. Linear systems are homogenous and additive. One more condition it satisifies is that if a signal is a sine or cos wave then the output only changes the phase or amplitude. Why? Why cant the output be a totally different output when a sine wave is given as an input?


texts: Python resampling and filtering question

texts:
Intro

Hi I'm looking for someone with some experience to walk me through some questions I have regarding resampling (downsampling) in Python.  Here is my problem:
I have accelerometer data that comes in irregularly at 50Hz.  Sometimes I get 50 samples a second, sometimes 1 every 3 seconds.  I also get velocity data that comes in at 1Hz also fairly irregularly.  My question is, how best to filter and align the data?

Problem

I have been using the pd.resample("1s").mean()  and pd.resample("1s").last() methods to downsample/decimate my data.  When using the .last() method I have been using a low pass filter with a cutoff of 0.5Hz prior to decimation.  As far as I can tell, the .mean() acts as it's own low pass filter, but I'm not sure if this is correct as there should be a difference between using a moving average filter and then resampling vs resampling into 1 second bins and then taking the mean.  I have seen in a lot of online dicussion people simple use pd.resample().mean() and that seems wrong to me.

Anyways, anyone with any insight would be greatly appreciated.
Cheers!


texts:
I am using IPP for FIR filtering over 1D data. I want to use the same input and filter it using two band pass filters. This is how I do it:

output1 = doFIRBandPassFilter(input,low1,high1);
output2 = doFIRBandPassFilter(input,low2,high2);

I want to combine those outputs so that the end result will show BOTH (low1,high1) and (low2,high2) filtered data.

If output1 and output2 are arrays - should I add them? multiply them? convolve them?


texts: Sound synthesis effects for natural sounds

texts:
I´m currently working on an accordion synthesizer and I would like to make it sound as natural as possible.

I really like how the following wave sounds:

https://dl.dropbox.com/u/20437903/onda%20acordeon.wav

The wave looks like this:

Looking at the wave it looks like it had harmonics added and substracted. Am I wrong?

I´m a bit new to sound synthesis and would be very grateful if anyone could please help me identifying the effects applied to the wave and some guidance to apply them. Thanks.


texts: Cooley-Tukey FFT: You don't have to zeropad to a power of 2?

texts:
Someone wrote "The algorithm that Cooley and Tukey presented in their classic paper (Math. Comp. 19 (1965), 297-301. http://dx.doi.org/10.1090/S0025-5718-1965-0178586-1) can be applied to any composite length.  The performance advantages are greatest for highly composite lengths, of which powers-of-2 are one example, and lengths of powers-of-2 result in other advantages on binary computers, so it has become a common misconception that the algorithm is only applicable to signals whose length is a power of 2."

Does that mean that when you DO use the Cooley-Tukey FFT You don't have to zeropad to a power of 2?
Take for example an image of 1920x1080. So, if you want to use the Cooley-Tukey FFT, you don't need to zeropad that to 2048*2048?


texts:
The research paper "Multidirectional Scratch Detection and Restoration in Digitized Old Images" says that,

  4.1. Preprocessing. The preprocessing step aims to enhance image features along a set of chosen directions. First, image is grey-scaled
  and filtered with a sharpening filter (we subtract from the image its
  local-mean filtered version), thus eliminating the DC component.

Now, I am Googling the term "Local Mean Filter" but there is nothing available like that.

Can anyone please provide me any reference of "Local-Mean Filter"?

What filter are they using for Sharpening?


texts:
I'm doing a project involving FMCW linear chirp signals.

Given the phase

$\phi = 2\pi\frac{k}{2}t^{2} + 2\pi f_{c} t + \theta$

where $k$ is the frequency increase rate.

When transmitting frequency modulated pulses continiously will the phase increase in this manner?

Or simply repeat in this manner?

What is the advantages and disadvantages of either one?


texts: How to separate the upwards-propagating from the downwards-propagating waves?

texts:
I've got a real 2D image of propagating waves in both directions. How can I separate this image into one with the upward-propagating waves and one with the downward-propagating ones?

This example has one wave in each direction. The x-scale shows time, the y-scale shows position, and intensity shows the amplitude of the wave at that position and time. The first image is my input.

However, I'm not very used to 2D signal processing, so the way I'm doing this is by just setting two quadrants of the FFT spectrum to zero and then reverse-transforming back. I got the other two images using the following algorithm:

X = fftshift(fft2( x )); % x is the input (image above)
X( 1:floor(size(X,1)/2),    1:floor(size(X,2)/2)   ) = 0; % null 1st quadrant
X( ceil(size(X,1)/2)+1:end, ceil(size(X,2)/2)+1:end) = 0; % null 3rd quadrant
y = ifft2( ifftshift(X) );
y = real( y ); % y is the output (image above)

It feels like such a barbaric method would introduce all kinds of artefacts, so I'm looking for feedback from someone more experienced.

Is there a better way to achieve this?
How should I "soften" the edges of the 0/1 "mask" I'm multiplying the FFT with?
The output is complex, so I'm discarding the imaginary part. Is that correct?
I'm also a bit worried about the fact that the noise will not be white anymore, but I suppose there isn't much to do about that except for low-pass filtering the whole image.


texts: Estimate sampling frequency given an array of values

texts:
Given an array x of size 2000, is there any way to calculate (or estimate) the sampling frequency using Matlab? The values are the only information I have about the signal. I know the signal has been sampled with Fs = 360 Hz but I want to know if there is any way to extract this information from these plots.

This is the signal plot, it is an ECG:

And this is the FFT:


texts: Inferring space domain signal from 2D DFT

texts: By just looking at the 2D Fourier Transform of a signal, can it every be known precisely which values in the space domain are zero?


texts: what is the use of signal processing in CT scans?

texts:
I'd like to ask about the use of signal processing in CT scans ?
would be great if you suggest Book names or courses about it !
 thanks in advance !


texts: Can I use the word amplitude instead of magnitude when I describe FFT bins? I dont see any similar word in my language.


texts:
What effects does have a different input such as impulse, step, sine sweep etc. on performing a system identification?

I mean what is up/down side, and which one is best, and why chose one rather than the other.

For now I see the sine sweep as the most reliable for determining a system since the other ones doesn't seem to hit nearly all modes of a system.

But what do you say? I am not an expert at all and I am a bit confused here.


texts: How $\Delta f=\frac{1}{2T}$ satifies the orthogonality condition?

texts:
Let
$$s_{ml}(t)=\sqrt{\frac{2E}{T}}\exp(j2\pi\Delta fmt)$$
where $T$ is the time-period of signal
$\Delta f$ is the frequency spacing
The text says that two signals $s_{ml}(t)$ and $s_{nl}(t)$ are orthogonal to each other for $\Delta f=\frac{1}{2T}$

Orthogonality means that $$<s_{ml}(t),s_{nl}(t)>=0, \text{when}~m\neq n$$
$<>$ is inner product

On further simplifying :
$$\int_{0}^{T}\frac{2E}{T}\exp(j2\pi\Delta f(m-n)t)~\text{d}t$$
which further simplifies to
$$\frac{2E \sin(\pi T \Delta f (m-n))}{\pi T \Delta f (m-n)} \exp(j\pi\Delta f(m-n)T)$$
which is precicely
$$2E~\text{sinc}(T \Delta f (m-n))\exp(j\pi\Delta f(m-n)T)$$
As per my understanding the aforementioned term should be zero for $\Delta f=1/T$ and not $\Delta f=1/2T$.
Where I am wrong?


texts: Converting mel spectrogram to spectrogram

texts:
I have a set of songs for which I extracted the STFT (Short-Time Fourier Transform) and used the magnitude spectrum $|S|$ to calculate the mel spectrogram by using a mel filterbank matrix $M$, so $X=\log(M\times |S|)$. I want to know is there any method to reverse this process, i.e. convert from the mel spectrogram back to the spectrogram. I performed some dimensionality reduction on the mel spectrogram, and reconstructed the mel spectrogram from lower dimensions. Now I want to regenerate the audio signal from the reconstructed mel spectrogram, so I guess first reconstruct the spectrogram and then the audio signal.

The problem is that the mel filter bank matrix is not a square matrix, since we the reduce the no of frequency bins, so inverse of $M$ cant be used like this : $ \hat{S}=M^{-1}\exp(X)$. So is there any way to generate the inverse mapping, like some inverse transfer function that can convert from $X$ to $S$?


texts: FMCW Radar: Need for Second stage FFT for velocity determination?

texts:
I have understood that for the purpose of Range determination, the samples of the beat frequency are first assembled per chirp. The beat frequency is directly proportional to the range of the reflecting object.

So if there are N chirps, there will in a 2D matrix be N Rows of M samples each. When the FFT is performed along each row (Chirp), the spectrum analysis yields the value of the beat frequency of the corresponding chirp. Let us assume that the beat frequency is represented by Bin B.

If the object is stationary, then all of the rows will have a significant value and nearly the same value in Bin B [Magnitude and Phase].

If the object is slowly moving away and assuming that each bin represents a large range, then Bin B continues to remain the dominant bin with possibly reducing magnitude and increasing phase.

Question hence is: Why perform a second round of FFT along a column for determining velocity of a range bin?

Why couldn't the individual phase of the FFT samples [Output of first round of FFT] along a column be taken and the velocity deduced from phase difference among bins of the column of a range bin?

After all, the phase of a received signal is proportional to the range of the object.


texts:
I'm writing a biquad filter program in C++, using the LBJ Cookbook as the source for the formulas. https://webaudio.github.io/Audio-EQ-Cookbook/audio-eq-cookbook.html
The lowpass filter that I've made seems to half work. The frequencies above the cutoff are being attenuated, however, it is causing aliasing similar to what happens with frequencies above Nyquist's limit.
double output, a0, a1, a2, b0, b1, b2, w, alpha;
double pi = 3.14159265359;
w = (2 * pi) * (cutoff / sampleRate);
alpha = (sin(w)) / 2 * q;
a0 = 1 + alpha;
a1 = -2 * cos(w);
a2 = 1 - alpha;
b0 = (1 - cos(w)) / 2;
b1 = 1 - cos(w);
b2 = (1 - cos(w)) / 2;

output = (b0 / a0) * inSamples[0] + (b1 / a0) * inSamples[1] + (b2 / a0) * inSamples[2]
  - (a1 / a0) * outSamples[0] - (a2 / a0) * outSamples[1];

I'm not sure where my error is, so any help is greatly appreciated!


texts: Why Does ASK Modulation Create Fourier Sidebands?

texts:
I know why analog amplitude modulation has side bands, it is related to (fc+fd) and (fc-fd). But what about DAM?
ASK(DAM) is a type of digital modulation, and there are only two state: carrier signal for "1" or nothing for "0". The expected behavior on a spectrum is seeing only one peak on the fc, but there are many side bands at fc+3d, fc-3fd, fc+5fd, fc-5fd etc.
Can someone explain that why we see these Fourier series pattern side bands?
fc: carrier frequency
fd: data signal frequency


texts:
I have a pressure transducer signal from a reacting mixture available here, which has a sampling frequency of 50kHZ. I have also obtained the spectrogram of the signal with NFFT=256, noverlap=4 as shown below:

I expected to see one high power frequency after pressure rise, instead of multiple frequencies (after t=0.03s). Why are multiple frequencies detected? Am I using the correct NFFT and noverlap parameters? How are NFFT and noverlap selected for this signal?


texts: Why is the Fourier Series a special case of the Fourier Transform and not the other way around?

texts:
I was reading a text book on the frequency domains in signal processing and my understanding is that the Fourier Transform considers signals that are a-periodic in time while the Fourier Series considers periodic signals in time (such as a sum of sinusoids).
But the book clearly mentions that the FS can be seen as a special case of the FT and that doesn't make sense to me.
My understanding is that the FT is a special case of FS where periodicity is just infinity. But the book says its the other way around without providing its own reasoning.
It makes no sense and I can't find answers online


texts: Non-uniqueness of continuous time sinusoids?

texts:
I have read the DSP literature and I get why:
$$e^{j(\omega+2 \pi m)n}=e^{j\omega n}$$

But can't we prove the same in continuous domain also in the following way:
$$e^{j\omega t}$$
$$e^{j(\omega+2 \pi m)t}$$
$$\left( e^{j\omega} \cdot e^{2j\pi m}\right) e^t$$
$$e^{j\omega t}$$

If my above reasoning is correct, then does not it mean that analog frequencies separated by $2\pi$ are same? But I know this is not the case.

Please help me in this case.


texts: Avoid saturated pixels in cross-correlation computation

texts:
I am doing template maching using Normalized Cross-Correlation for measuring similarity between template and image.

However, some images from my dataset contain saturated pixels (overshoot) which should be avoided in the measurement.

How to update the cross-correlation computation to avoid saturated pixels?


texts: How to measure SINAD of an audio signal?

texts:
I tried to measure the SINAD and SNR of a recorded audio signal using MATLAB. The following is the process I did to measure the 2 parameters:
Firstly, I generated an FM modulated signal(1 kHz audio tone) at a particular frequency using a signal generator and connected it to an RTL-SDR. The RTL-SDR is connected to a PC with a loaded SDR# software (even GNURadio can do the operation). After tuning the RTL-SDR to generated RF signal, I was able to receive the 1 kHz audio tone. I recorded the sample stream for a small duration (1 minute) and loaded the recorded data in MATLAB. I used the in-built SINAD and SNR function to calculate the sensitivity and signal to noise ratio. Surprisingly, I receive the same values for both SINAD and SNR. But, When I look at the graphs for SINAD as well as SNR, it appears to be that both parameters are calculated wrongly.

From the graphs, It appears that both SNR and SINAD should have higher values like more than 60 dB.
If you have any ideas, on how to calculate these parameters in MATLAB, please let me know.


texts:
I have just started my first proper course in DSP, and am a little confused about a statement in my textbook.  The book states:

"Notice that when the Nyqust-frequency is defined as $\pi/\Delta t$, we will with $\omega_m = 2 \pi m/N \Delta t$ where $m = 0, \pm 1, \pm2, ...$ have $\omega_m \neq \pi/ \Delta t$, unless $N$ is even.  If we are working with symmetric functions, it can be wise to use an odd number of sampling points."

I see the logic in the first sentence here, but don't understand how it follows from this that it is wise to use an odd number of sampling points for a symmetric function. This statement is not elaborated upon further. Can anyone explain this logic to me?  Why it is wise to use an odd number of sampling points for a symmetric function?  If anyone could explain this, I would greatly appreciate it!


texts:
Is it possible to differentiate between FIR and IIR Filter just by knowing the poles and zeros of the system?

If yes, how?
If not, why?

Thank you.


texts: What is the tightest known bound on the reconstruction error in compressed sensing?

texts:
I am specifically thinking about the reconstruction error of L1-minimization from compressed measurements with noise.
I know a bound from (8) in The restricted isometry property and its implications for compressed sensing, but I was wondering if there is anything newer and or tighter that I am not aware of. Additionally, I would like to relate $\delta_{2s}$ to the size of e.g. a Gaussian measurement matrix. It seems I can do that from Lemma 3.1 and (3.22) in Decoding by linear programming, but how tight is (3.22) and can I get closer?


texts:
I just can't get my head around Fourier transform and convolution in 2D. I am trying to implement image convolution using fast Fourier transform (in julia). So the first thing I need to do is to pad the image with zeros. Let's assume that the image is $N\times P$ and the kernel $s \times q$, so I would pad the image so that it is of size $(N+s-1)\times (P+q-1)$.
(N, P) = size(image)
(s, q) = size(kernel)
image_padded = zeros(N+s-1, P+q-1)
image_padded[1:N, 1:P] = image

Then I need to make the kernel the same size as the image, so I pad it with zeros, too, but this time, I have to put the kernel in the center of the huge block of zeros.
kernel_padded = zeros(N+s-1, P+s-1)
Npul = Int((N+s-2)/2) #both image and kernel have odd sizes
Ppul = Int((P+s-2)/2)
spul = Int((s-1)/2)
qpul = Int((q-1)/2)
kernel_padded[Npul-spul+1:Npul+spul+1, Ppul-spul+1:Ppul+spul+1] = k

Then I take fft of both padded image and padded kernel, take their Hadamard product and use inverse fft to get the result back to image domain.
real.(ifft(fft(image_padded).*fft(kernel_padded)))

But the result does not look as I would expect, there is like a cross in the middle of my resulting image. Actually, it looks like there is one quater of the image in each corner but rotated so that the border effect creates the dark cross in the middle. The only way I can get rid of it is to pad the image and kernel to $(N+N+1)\times(P+P+1)$.
So do I understand the process wrong or is something wrong with my code? Thanks for any advice.


texts: Should the transmitted signal be continuous in FMCW?

texts:
I am a novice in this field and working on respiration detection through near ultra-sound between $18 \textrm{kHz}$ and $20 \textrm{kHz}$.

In the referece article and other radar use articles, the generated signal was transmitted continuoulsy such as the picture below:

My question is what if the interval would be allowed between the signals such as the picture below:

I would like to know why the signla be tranmitted continuously and whether it is possible to be transmitted with intervals.

Thank you in advance.


texts:
I know that the error probability for binary antipodal signaling is $Q(\sqrt{2E_b/N_0})$ with $Q$ the tail function.

I was wondering if when used in a m-ary scheme how the probability error would change. What I mean by that would be a signaling where we still have only 2 antipodal codewords, but in a ${\rm I\!R^m}$ dimension domain, i.e. $c_1$ = +$\sqrt{E_b/m}$(1,..1)
$c_2$ = -$\sqrt{E_b/m}$(1,..1)

Thank you for reading, and have a great day.


texts: A simpler method or more descriptive answer to the Fourier Transform

texts:
I'm trying to do a CT Fourier Transform of these two signals
$$x_1(t)=e^{−a(t−1)} \cdot u(t−1)$$
and
$$x_2(t)=e^{−a(t−1)} \cdot u(t)$$
Where $a$ is any real number, and $u(t)$ is the unit step function and that is multiplication, not convolution.

Picture of my work and answers here.

The answers listed should not contain the extra
$$t \cdot u(t-1)$$ and $$t \cdot u(t)$$

Okay, thanks to Dilip i now have my answers to look like
$$ x_1(F)=\frac{e^{2a+j2\pi\cdot F}}{a+j2\pi\cdot F}$$
and
$$x_2(F) =\frac{e^a}{a+j2\pi\cdot F} $$
That does seem more elegant than what i previously had, I feel like I could possibly manipulate $x_1(F)$ into a sinc function, but I don't think it's necessary.


texts: Why do wavelets integrate to 0 and how do they match a signal?

texts:
I have been reading about the Wavelet transform recently and its relationship to the Fourier transform.  From what I understand the wavelet transform represents signal data with many short-lived functions (wavelets).  The major difference between the 2 is that the Fourier does not sample in time, but only in frequency.
 I have  couple questions about the wavelet method:

1) Why is it required that wavelet functions integrate to 0?

From my understanding, wavelets are shifted along the time domain and compressed or stretched to match the frequency of the signal at that point in time.  So wavelets are able to capture frequency at different places in time.

2) How does the wavelet capture peaks that are higher or lower in the original domain?  For example see the image below that I have simulated in R.  The dashed line represents the wavelet transformed approximation.  How  is this dashed line able to capture peaks that are higher or lower if it is only shifted and compressed or stretched in the x direction?


texts: Estimate gap between two adjacent images

texts:
I have two images of a panoramic view but there is a small gap (<24 pixels) in the horizontal dimension between them. I would like to do interpolation/inpainting to fill the gap, but the exact size of the gap is unknown.

Is there a way to estimate the best value for the gap? My first thought is to detect lines/curves in the images (Canny edge detection et al) and somehow try to connect the lines at the adjacent edges, and see which value for the gap gives the smoothest continuation. However I don't known where to begin, or if there is a better/simpler solution.


texts: Predicting Positional Displacement with a Bivariate Gaussian Distribution

texts:
In my computer-vision problem, a blob can move either forward or backwards (with a slight chance of moving sideways) along its directional axis. I am seeking a bivariate model because I need to determine whether a blob's new position (with respect to its old one) is "feasible," or "likely." However, the blob can not move extreme distances very quickly.

Thus, there is a high probability that the blob will not move very far if it does move (but still a sufficient amount such that the probability distribution will be elongated along one axis [not the orthogonal]).

I am searching for a model that has high probability values within a smaller region, but has a larger encompassing region with small (and distinctly different) probability values.

Thank you for your help!


texts: Why does increasing FFT length (narrowing bandwidth) not decrease Noise Power per bin?

texts:
I was under the impression that zero-mean white guassian noise had a constant Power Spectral Density which means that smaller bandwidths should decrease the power in that band. I'm getting surprising results for average power in each FFT bin, it's staying the same irrespective of how long I make the FFT. Am I not measuring power correctly?

Here is a reduced for of my MATLAB script (should be Octave-compatible). Please excuse my odd habit of time vectors going up-down and frequency and misc ones going left-right.

%pkg load communications
% ^UNCOMMENT^ if using OCTAVE
close all

Q=8;        % Lowest power of 2 to iterate though in FFT loop
R=12;       % Highest power of 2 to iterate though in FFT loop
N=1;        % Linear PSD of noise
c=0;        % Loop counter

K=2.^(Q:R)  % Loop vector

for k=K
    n=sqrt(N).*(randn(k,1)+j.*randn(k,1))./sqrt(2); % complex noise
    P=(abs(fft(n)')).^2./k; % Noise Power in frequency domain
    c=c+1;
    e(c)=mean(P); % Mean power measurement across bins
    %v(c)=var(P);
    figure
    plot(10.*log10(P))
end
e
%v

I feel like I'm missing something fundamental about this. Does FFT binning not truly band-limit signal power because of the way the algorithm works? I'm basically trying to get an output analogous to a spectrum analyzer in that is shows something proportional to true power.


texts: confused in practical example of shift invariance?

texts:
I am reading book of w smith ,"The scientist and engineer guide to digital signal processing" 2 edition
I am trying to understand topic of shift invariance from theory given on pg 92

I have attached screen shot and i have indicated the example with in red marks.

I am confused in whether shift in variance ,author is talking about is in regard to which system?telephone line or digital filter?
What is our target system in this example that needs to be shift invariant?Is it telephone line or digital filter?


texts: How to implement a gradient based Hough transform

texts:
I am trying to use the Hough transform for edge detection, and would like to use gradient images as the basis.

What I have done so far, given the image I of size [M,N] and its partial derivatives gx, gy, is to calculate the gradient angle in each pixel as thetas = atan(gy(x,y) ./ gx. Similarly I calculate the gradient magnitude as magnitudes = sqrt(gx.^2+gy.^2).

To build the Hough transform, I use the following MATLAB code:

max_rho = ceil(sqrt(M^2 + N^2));
hough = zeros(2*max_rho, 101);
for x=1:M
    for y=1:N
        theta = thetas(x,y);
        rho = x*cos(theta) + y*sin(theta);

        rho_idx = round(rho)+max_rho;
        theta_idx = floor((theta + pi/2) / pi * 100) + 1;
        hough(rho_idx, theta_idx) = hough(rho_idx, theta_idx) + magnitudes(x,y);
    end
end

The resulting Hough transform looks plausible (see http://i.stack.imgur.com/hC9mP.png), but when I try to use its maxima as edge parameters in the original image, the results look more or less random. Did I do something wrong in constructing the Hough transform?

UPDATE: I had a stupid mistake in my code: rho was calculated as x*cos(theta)+y*cos(theta) instead of x*cos(theta)+y*sin(theta). That is, I was using two cosines instead of a cosine and a sine. I have edited the code above and the new resulting image is below. This did not give much better edges though.

@endolith: To plot an edge, given a maximal value in the hough-matrix at rho_idx, theta_idx, I translate the indices to rho,theta values:

theta = (theta_idx -1) / 100 * pi - pi / 2;
rho = rho_idx - max_rho;

Finally I plot the edge as y= (rho - x*cos(theta)) / sin(theta).


texts: use wavelet transform to extract frequencies from given signal


texts: How is the topo / scalp map from an EEG Signal with ICA created?

texts:
I am lately experimenting with an EEG and found this Library: MNE Python.
I am using an ICA to find and correct artifacts - this is working fine.

It's now possible to plot a topo / scalp map for the ICA (EEGLab and other also can do this).

Now I wonder how this is calculated?

I assume it takes the topographic information of the EEG sensors and combines them with the ICs for each channel. Is this correct?
Are there any sources explaining the process in general?

Thanks in advance


texts: Software implementation recipe of cascaded Biquad IIR filter

texts:
I have been given the task of measuring the performance of Biquad IIR filters on a new CPU that we are evaluating.

What I know is that the Biquad filters can be cascaded to form higher order filters.

I drew the following based on some light reading to form a fifth order IIR filter.

For every new sample that arrives, do I have to simply execute the 3 equations for each of the stages sequentially

y1 = b0*x + ..;
d11 = ..;
d12 = ...;
y2  = ...;
d21 = ..;
...
...
d31 = b1*y2...;

The 5 coefficients for each of the stages are available and so are the initial differences (Dx1 and Dx2). So should I just sequentially execute the stages for each sample? Is that all to Biquad IIR software implementation? I would be grateful if you can tell me what else is in the recipe.


texts: Implementation of a filter in a software application

texts:
I have to embed a bandpass filter into a software application.

All the signal processing will be done in 32 or 64-bit double precision.

The filter coefficients for an IIR Butterworth filter are computed offline in Matlab.

What is the best way to implement the filter in the application?

Right now my plan is to implement the filter as cascaded second order sections in transposed direct form II, but I have no practical experience (and literature on the practical part of implementation is also kind of sparse), so I'm just kind of guessing.

For example: Second order sections are used to avoid numerical problems, but do I really need to take care of that when the processing is done in floating point precision, also which form should i chose for the realization of the filter?

I would be thankful for any tips or literature hints.


texts:
I have an image which has a vertical line which is 2 pixels wide and 4 pixels long and the remaining image is black.

Now I want to find the edges in the vertical direction and for that I use the sobel operator

sobel = { {-1, 0, 1 }, { -2, 0, 2 }, { -1, 0, 1 } }

and get the following thresholded image

but when I use 180 degree rotated version of the above operator such as
sobel = { {1, 0, -1 }, { 2, 0, -2 }, { 1, 0, -1 } }

I get the following image

Now my question is

Using a 3 * 3 operator, do we get 2 pixels wide edge as in the previous two images? Shouldn't we get 1 pixel wide edge? and what about the edge width in case of 5 * 5 or 7 * 7 operators ?
What is the difference between the above said operators i.e normal sobel operator and its 180 degree rotated version. Where is each one used? In this case which of the above should be used?


texts: Typical parameters and stability of Chebyshev filters

texts:
I design Chebyshev first and second type filters using scipy.signal.iirfilter.

What are typical values of maximum ripple in the passband and minimum attenuation in the stopband? I struggle to find optimal values of those parameters.
I know that IIR filters are often unstable, but is it possible to make it stable by changing values of maximum ripple in the passband and minimum attenuation in the stopband? If not, are there any ways to make it more stable, say, for the 6th order filter?


texts: How to separate a low amplitude signal of same frequency but different phase

texts:
I have a signal consisting of two sine waves of same frequency but have different phase (and also a lot of HF noise). Their amplitudes differ by about $120\textrm{ dB}$. How to I go about separating them? SNR is about $130\textrm{ db}$ but its frequency is far higher than my signal, i.e. $f(\textrm{noise}) \gg  f(\textrm{signal})$.

Any special technique? Also I am really just a beginner, so extra points for those who give a detailed explanation.


texts:
I have an interlaced image and I like to make a deinterlaced image out of it. Currently I just take average of two rows to create the missing row. Is there a better way to do this? What is the state of the art?

I don't have the previous or the next frames, just this one.


texts:
I have k-space data related to an image. k-space data is made of Kx, Ky, signal arrays of 4100x1. What I have done for image reconstruction is:

$$I_{m,n}=\sum_i Signal_i \space e^{2\pi j(\frac{K_{x}^i \space m}{4100} + \frac{K_{y}^i \space n}{4100}) } $$

However, I cannot reconstruct the image correctly. What do you think is the mistake? Thanks in advance!


texts: Why does the general formula for DFT of real sinusoids not apply here?

texts:
In this course (Coursera: Audio Signal Processing for Music Applications) the professor derives a general equation for finding the DFT for real sinusoids:

  DFT of real sinusoids
  \begin{align}
x_3[n]&= A_0\cos\left(2\pi k_0 n/N\right)=\boxed{\frac{A_0}{2}e^{j2\pi k_0 n/N} +\frac{A_0}{2}e^{-j2\pi k_0 n/N}}\\
X_3[k]&=\sum_{n=-N/2}^{N/2-1} x_3[n]e^{-j2\pi k n/N}\\
&=\sum_{n=-N/2}^{N/2-1}\left(\frac{A_0}{2}e^{j2\pi k_0 n/N} +\frac{A_0}{2}e^{-j2\pi k_0 n/N}\right)e^{-j2\pi k n/N}\\
&=\sum_{n=-N/2}^{N/2-1}\frac{A_0}{2}e^{j2\pi k_0 n/N}e^{-j2\pi k n/N} +\sum_{n=-N/2}^{N/2-1}\frac{A_0}{2}e^{-j2\pi k_0 n/N}e^{-j2\pi k n/N}\\
&=\sum_{n=-N/2}^{N/2-1}\frac{A_0}{2}e^{-j2\pi\left(k-k_0\right) n/N}+\sum_{n=-N/2}^{N/2-1}\frac{A_0}{2}e^{-j2\pi\left(k+k_0\right) n/N}\\
&=N\frac{A_0}{2}\text{ for }k=k_0,-k_0; 0\text{ for the rest of }k
\end{align}

Keeping this in mind, given a discrete signal $x[n]$ that has 4 samples given as: $[1, -1, 1, -1]$

the equation from the slide can be used to represent our signal as:
\begin{align}
x[n] &= {A_0\over2}e^{j2 \pi k_0 n/N} + {A_0\over2}e^{-j2 \pi k_0 n/N}\\
&= {e^{j 4 \pi  n/4}\over2} +  {e^{-j 4 \pi  n/4}\over2}\\
&= {e^{j \pi  n}\over2} +  {e^{-j \pi  n}\over2}
\end{align}

And now again, calculating $X[2]$ with the equations from the same slide:

\begin{align}
X[2]&=\sum_{n=-N/2}^{N/2-1} \left({e^{j \pi  n}\over2} +  {e^{-j \pi  n}\over2} \right)  e^{-j2\pi kn/N}\\
&=\sum_{n=-N/2}^{N/2-1} \left({e^{j \pi  n}\over2} +  {e^{-j \pi  n}\over2} \right)  e^{-j4\pi n/4}\\
&=\sum_{n=-N/2}^{N/2-1} \left({e^{j \pi  n}\over2} +  {e^{-j \pi  n}\over2} \right)  e^{-j\pi n}\\
&=\sum_{n=-N/2}^{N/2-1} {e^{j (\pi n - \pi n)}\over2} +  {e^{-j 2 \pi  n}\over2}\\
&= \frac N2 + \frac N2\\
&= N
\end{align}

But this contradicts the last line on the slide that says $X[2]$ should be $N{A_0\over2}$.

Are the slides wrong?
I'm also having trouble figuring proving/understanding the following equality from the lecture slide:

  $$\sum_{n=-N/2}^{N/2-1} {A_0\over2} e^{-j2\pi(k-k_0)n/N} + \sum_{n=-N/2}^{N/2-1} {A_0\over2} e^{-j2\pi(k+k_0)n/N} = N\frac{A_0}{2}\text{ for } k=k_0, -k_0$$

PS: From the same course, some slide from a another lecture show these equations that are consistent with my derivation (see equation for $<x,s_2>$)

  DFT: scalar product:
  $$
\boxed{<x, s_k>=\sum_{n=0}^{N-1}x[n]s_k^{*}[n]=\sum_{n=0}^{N-1}x[n]e^{-j2 \pi k n/N}}
$$
  Example:

  $$x[n]=[1, -1, 1, -1]; N=4$$

  \begin{align}
<x, s_0>&=1*1 +(-1)*1 + 1*1 + (-1)*1 =0\\
<x, s_1>&=1*1 +(-1)*(-j) + 1*(-1) + (-1)*j =0\\
<x, s_2>&=1*1 +(-1)*(-1) + 1*1 + (-1)*(-1) =4\\
<x, s_3>&=1*1 +(-1)*j + 1*(-1) + (-1)*(-j) =0\\
\end{align}


texts: Spectrum analysis : how to represent FFT results

texts:
I'm trying to program a spectrum analyser in C++ but I'm faced with a conceptual problem : how to map the results of the FFT to an image ?

Let me detail :
I'm doing all this in a dynamic way.

I have chunks of audio data as float [-1;1] arriving to my class.
I put them in a circular buffer which has the same size as my FFT window
I apply a Hann(ing) windowing function
I apply an FFT to the result.

This part works great.

Now the result of my FFT is a two dimensional array (my library automatically does the sqrt(re^2 + im^2) part) with bin indexes and amplitudes.

My goal is to draw this result from time to time into a window. The problem is, the window has way less pixels than there are bins (if my buffer size is 65536, I have 32768 bins), so I have to do some kind of "averaging".
I could do a mean, a root square mean, take the maximum, etc etc ... I have no idea what is the best method to use ? The application of this will be music, more specifically, detect frequency masking by plotting several semi-transparent sepectrums on the top of each other.

Any advice on a good way to do that would be much appreciated


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


texts: Understanding BFSK Modulation (MATLAB/Simulink)

texts:
Could someone explain to me BFSK Modulation on the Simulink model below?

The simulation above simulates a binary FSK system with a correlation (or matched filter) receiver.

I have to understand the general idea about it.

Could someone explain this?
Moreover I cannot understand the meaning of correlator and tone frequency fn the above simulator?


texts: Is there an analytic expression by which it is possible to calculate overshoot of step-response of second order IIR digital filter? I have tried to search through Mitra's signal processing book and through some other materials, but wasn't able to find it.


texts: MRI image reconstruction - reduced scan pecentage

texts:
I have an MRI kspace from a siemens 1.5 T scanner with reduced scan percentage ( sampled truncation by not acquiring the most peripheral lines in kspace). The kspace plot in matlab look like below.

I can do zero filling in this kspace to the required resolution and reconstruct images.

Other than zero filling, is there another method which can estimate missing lines and reconstruction images.


texts: LMS algorithm for modeling step-size ambiguity!

texts:
Behrouz Farhang-Boroujeny in his adoptive filters, 2nd ed., p. 155 ,told:

  It is sufficient for stability: 1/(3*tr(R)).

But his book have attached mfile for modeling:

%   Modeling
% Last updated on April 28, 1998
%
itn=input('\n No. of iterations?      ');
sigman2=input('\n Variance of the plant noise?      ');
sigman=sqrt(sigman2);
wo=input('\n Plant impulse response (vector, w_o)?      ');
a=size(wo);
if a(1)<a(2)
   wo=wo';
end

N=input('\n Length of the model (N)?      ');

h=input('\n Coloring filter impulse response (vector, h)?     ');
a=size(h);
if a(1)<a(2)
   h=h';
end

Misad=input('\n Misadjustment (e.g., 0.1 for 10%) ?     ');
mu=Misad/(N*(h'*h));

You can see:

mu=Misad/(N*(h'*h))

And it is not matching to his formula in his book ...
I can't understand where the hell it is coming from?

More information about problem:

Scheme of problem:(Note the source of image is irrelevant to my question, it's only image)

Core of algorithm: (Note, it is average over runs for smooth learning curve)

for k=1:runs
    x=filter(h,1,randn(itn,1));% the randn creates input sequence that colorize by filtering through filter with h impulse response
    d=filter(wo,1,x)+sigman*randn(itn,1); % Create desired by sum of noise and y(main plant output)
    w=zeros(N,1);

    for n=N:itn
        xtdl=x(n:-1:n-N+1);
        e=d(n)-w'*xtdl;
        w=w+2*mu*e*xtdl;
        xi(n)=xi(n)+e^2; % Creation of learning curve, gradient(e^2)=-2e(n)x(n)
    end
end
xi=xi/runs;

Note: 1/(3*tr(R)) make algorithm unstable ...!!! But without colorizing filter this will converge.


texts:
Suppose I have two signals $f,g$ which are orthogonal in the sense of the usual inner product on $\mathbb R^d$. Then if I'm not mistaken these signals have disjoint spectral support, i.e their fourier transforms have disjoint support. Suppose further that I know the spectral supports aswell. How can I make use of this information in order to sample only one of the signals, i.e, sample so I don't "see" the other at all?

The trivial solution seems involve filtering out the irrelevant frequencies using an ideal filter, but that's not an option; I am restricted to "physical" solutions.

It's fine to assume $f,g$ are well-behaved.


texts:
I am working on a FMCW radar and I want to sample the I/Q signals from the 2 receiving antenna (total 4 channels) simultaneoulsy and send the digital data to computer to extract range and velocity information using FFT.
The below diagram is a overall plan.

In the diagram I have shown a ADC with 1MSPS but thats not based on any calculation, just a rough idea as the maximum frequency component is 500kHz. To make it simplier and test for short range I planned to limit the frequency to 100 kHz.

Now my question is how should I senlect the sampling rate? I don't think 2*fmax (Sampling Theorem) is enough for this purpose. The output from the radar frontend are I/Q signals, band limited from 30Hz to 100kHz. My target is to digitize the signal to the best for proper analysis. Should I consider the DAC (used for tuning the frontend) sample number calculated from the minimum and maximum voltage for tuning, to calculate my sampling rate?


texts: Time invariance from convolution integral

texts: I wanted to prove, time invariance property of system by convolution integral(i.e equal time shift in input result in equal shift in output) but as calculation shown in figure, I got unexpected result, so what I'm missing?


texts: How to express linear convolution using positive frequencies from channel and symbol DFT vectors?

texts:
Sorry for the long question, but I'm struggling to implement something related to this question. Your help would be appreciated.

Background

Notation: In the sequel, the uppercase letters represent frequency domain vectors, while lowercase represent time-domain.

Suppose I have two $N$-length hermitian symmetric DFT vectors $\mathbf{X}$ and $\mathbf{H}$, one representing a symbol and the other a channel through which the symbol goes. Thus, the received symbol in the time domain could be expressed as:

\begin{align}
\mathbf{y} = \mathbf{h} * \mathbf{x}
\end{align}
where $\mathbf{h}$ is the channel impulse response, $\mathbf{x}$ is the IDFT of $\mathbf{X}$ and "$*$" denotes the circular convolution.

This result could also be obtained by:

\begin{align}
\mathbf{y} = \tilde{H} \mathbf{x}
\end{align}
where $\tilde{H}$ is the circulant matrix.

Now, considering only the positive frequencies from both $\mathbf{X}$ and $\mathbf{H}$, that is, $\mathbf{X_+} = \{\mathbf{X(1), X(2), \cdots, X(N/2+1)}\}$ and $\mathbf{H_+} = \{\mathbf{H(1), H(2), \cdots, H(N/2+1)}\}$, it is possible to obtain $\mathbf{y_+}$:

\begin{align}
\mathbf{y_+} = \mathbf{h_+} * \mathbf{x_+}
\end{align}
where $\mathbf{x_+}$ is the IDFT of $\mathbf{X_+}$ and $\mathbf{h_+}$ the IDFT of $\mathbf{H_+}$, both having size $(N/2 + 1)$.

Then, $\mathbf{y}$ could be perfectly reconstruced from $\mathbf{y_+}$, since its IDFT:

\begin{align}
\mathbf{Y}_+ = \mathbf{H_+}\cdot\mathbf{X_+}
\end{align}

The $\mathbf{y_+}$ vector can also be obtained by making a circulant matrix with the IDFT of $\mathbf{H_+}$, which is a complex vector (because $\mathbf{H_+}$ does not have hermitian symmetry), and multiplying the new circulant matrix by the IDFT of $\mathbf{X_+}$, which is also complex.

Goal:

The ideal received symbol $\mathbf{y}$ is of my interest. However, I also need to know the intersymbol interference occurring when $\mathbf{x}$ is linearly convolved (not circularly) with the channel. I want to know this using $\mathbf{x_+}$ and $\mathbf{h_+}$, since I want to pre-equalize the symbols and using $(N/2 +1)$ vectors requires less operations.

Question

Using the IDFT of the positive frequencies from the original hermitian symmetric DFT vectors  as the actual transmit symbols, how could the result of a linear convolution be recovered perfectly?
The circular convolution is perfectly recovered.

I currently think it is not possible, because the interference will be different.

Sample code

The circular convolution is perfectly recovered:

N=8;

x = rand(N,1); % Time-domain random symbol
h = rand(N,1); % Time-domain random channel
X = fft(x); % Symbol FFT
H = fft(h); % Channel FFT
Y = H.*X; % Ideal Received Symbol (Ideal because of circular convolution)

% Time-domain complex vectors from positive frequencies
x_plus = ifft(X(1:(N/2 +1)));
h_plus = ifft(H(1:(N/2 +1)));
% Time-domain complex received symbol
y_plus = cconv(x_plus,h_plus,(N/2 +1));
Y_plus = fft(y_plus); % FFT of the complex received yields positive freqs.

% Comparison
[Y_plus; flipud(conj(Y_plus(2:end-1)))]
% Should be equal to:
Y

It doesn't work for linear convolution:

N=8;

x = rand(N,1); % Time-domain random symbol
h = rand(N,1); % Time-domain random channel
X = fft(x); % Symbol FFT
H = fft(h); % Channel FFT

y = conv(h,x); % Linear convolution
y = y(1:N); % Remove convolution tail
Y = fft(y); % Ideal Received Symbol (Ideal because of circular convolution)

% Time-domain complex vectors from positive frequencies
x_plus = ifft(X(1:(N/2 +1)));
h_plus = ifft(H(1:(N/2 +1)));
% Time-domain complex received symbol
y_plus = conv(x_plus,h_plus); % Linear convolution
y_plus = y_plus(1:(N/2+1)); % Remove convolution tail
Y_plus = fft(y_plus); % FFT of the complex received yields positive freqs.

% Comparison
[Y_plus; flipud(conj(Y_plus(2:end-1)))]
% Should be equal to:
Y


texts:
I am working on a project, for which I am finding the onset times of an audio and then extracting the MFCC coefficients for all the times.

Thus if I have n onset times, I would have n feature vectors, corresponding to the MFCC coefficients at those points.

Now I want to characterize a "bol" of the tabla by the the waveform in between the start of one onset and the start of the next onset. This would essentially mean that for one onset time, I would have to take the MFCCs feature vectors at multiple points, and I would have a matrix of onset values for each onset point, instead of just a single vector at each point. Also, the matrix wouldnot always be of the same size, since the gap between 2 onsets might not always be the same.

How do I combine these multiple MFCC vectors to one, or how to I use this matrix for machine learning?


texts: FFT followed by SVD leads to topography?

texts:
I am implementing the localizer method from this paper. One of the steps is not hard to understand and implement, but I don't understand why it is applied or what is the rationale behind it.
Each data epoch is a 60x500 array corresponding to 500 points on 60 channels sampled at 500 Hz. The paper states:

Fourier transform with 1 Hz resolution was applied to determine the
frequency in 8–13 Hz band that had
maximum power per Hz (frequency of maximum power (FMP)). Since the
amplitude-modulated sound occurred in 40 Hz, the Fourier transform was
applied to extract real and imaginary parts of 40 Hz as well as in FMP
for each epoch.

8-13 Hz corresponds to the alpha band. The amplitude-modulated audio sound occuring in 40 Hz is called ASSR. My python implementation of this step:
for i, epoch in enumerate(epochs):
    # epoch.shape = (60, 500)

    # Util function to return the fq associated to each bin
    fft_freq = np.fft.rfftfreq(epoch.shape[1], 1.0/fs) # 0 to 250, step 1
    # FFT
    fft_val = np.fft.rfft(epoch, axis=1) # (60, 251)

    # Keep the FFT for alpha band and ASSR
    alpha_band = np.where(np.logical_and(fft_freq>=8, fft_freq<=13))[0]
    alpha_band_fft = fft_val[:, alpha_band] # (60, 6)

    ASSR = np.where(fft_freq == 40)[0]
    ASSR_fft = fft_val[:, ASSR] # (60, 1)

    # Find FMP
    psd = np.square(np.abs(alpha_band_fft)) # (60, 6), equivalent to a periodogram
    FMP = 2 # FMP determine in 0,1,2,3,4,5 with 0 corresponding to 8 Hz and 5 to 13 Hz
    FMP_fft = alpha_band_fft[:, [FMP]] # (60, 1), In this example, I selected 10 Hz

As I am not yet sure how the FMP is chosen, I simply took 2 (corresponding to 10 Hz) as an example here.

We applied SVD on real and imaginary parts on FMP and 40 Hz separately
to calculate two different topographies of the brain in each epoch.
The first column of u for 40 Hz is the topography of
40 Hz and the first column of u for FMP is the
topography of FMP.

The relation between u and SVD is explicitly given in the paper, just in case. I used numpy for the implementation.
u, s, v = np.linalg.svd(ASSR_fft, full_matrices=False)
u, s, v = np.linalg.svd(FMP_fft, full_matrices=False)

With the parameter full_matrices set to True, the shapes returned are: u.shape = (60, 60), s.shape = (1, ),  v.shape = (1, 1). However, since only the first column from u is interesting, I have set full_matrices to False to reduce computation and to return an array u of shape (60, 1).
Question: What is the information extracted here? Why is the first column of u the topography of a given frequency? I do not get what the SVD is achieving here.


texts: Calculating 1/3 Octave Spectrum from FFT / DFT

texts:
I am not often on this forum and I am not an expert on the subject.

I struggle with the theory of FFT / DFT and the 1/3 octave spectrum.

Assume I have a DFT analysis of a given signal.
It (the DFT analysis) consists of many equidistant frequency bins that contain the corresponding amplitude (RMS or Peak).
I now want to calculate a 1/3 octave spectrum which has different frequency bins that are not equidistant.
Assume that one of the 1/3 octave bins reaches from 100 Hz to 200 Hz and the bin size of my DFT is 10 Hz.
Assume the amplitude of the DFT is already in RMS.

Here is what I do at the moment (in Matlab) but it seems that it's  not correct:

Adding all DFT amplitudes within 100 Hz and 200 Hz in linear scale (Pascal in my case).
After adding the RMS amplitudes (within one 1/3 octave bin) I just calculate corresponding the Decibel (dB) value - that's it.

Is there more to it? Can you point me in the correct direction?

  I do not "want" (I have to mimic another software) to use so-called
  1/3-Octave-Band Filter Banks like in this Matlab example.

Related

1/3 octave spectra from fft


texts: Have you ever seen logarithmic burn-in on a Kalman filter before?

texts:
I'm reviewing a paper that uses a multiplicative factor $\ln i$ where $i$ is the number of steps since the Kalman filter was initialized.  The idea is to slowly build confidence as the filter sees more data.

The authors replace the typical Kalman gain $K$ with one weighted by this multiplicative factor $K_{new}$:

  $K_{new} = ( \ln i) K$

Have you ever seen such an approach used before? I'd like to give the authors a reference to cite to tie them into the literature on this subject if it exists.


texts: Conversion from Floating point representation to Fixed point representation

texts:
I have written a matlab function for a constellation mapper and I want to create a Fixed-point version of that code for implementation on FPGA. I know matlab has a specific toolbox for fixed-point conversion but i want to do this manually. My specific question is:

1) which parts of the code would require this conversion? i mean do all inputs, outputs, temporary  variables need to be converted to a fixed point format ?

2) What are the steps that I should follow to do this conversion? i read that converting from floating point to fixed is no more, no less than scaling numbers. but I'm not quite sure how I go about doing this. Also, I want this fixed point representation to be of a custom length, not necessarily equal to standard lengths like 8, 16, 32-bits .

PS: I don't want to compare the performance of the fixed point code and floating point code in the time being


texts: Choosing the length of a signal when calculating the FFT

texts:
I'm trying to transforming an already know signal from time domain to frequency domain using fft function in MATLAB. And I got different values of amplitude when I change the value of length of signal $L$. I'm therefore wondering if there is a rule to choose the length of signal when calculating FFT.

I have an acoustic signal got from a sensor and need to analyze it (I have no idea about its frequency range, amplitude, and so on.) I want to transform it into frequency domain using fft function in MATLAB. And I don't know how can I choose the length of signal.

Could somebody please help?


texts: How to find H(z) and H(k) from a given causal function

texts:
Consider the causal function,
$y[k] = 2x[k] - 40x[k - 1] + 10y[k - 1]$  $16y[k - 2]$;
where $y[k]$ is the output and $x[k]$ is the input. Assume that the system is initially at rest.
Please someone help me find $H(z)$ and $H(k)$ from given equation. I have an assignment due to tomorrow